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CONVERSION  FACTORS,  NON-SI  TO  SI  (METRIC) 

UNITS  OF  MEASUREMENT 

Non-SI  units  of  measurement  used  in  this  report  can  be  converted  to  SI 
(metric)  units  as  follows: 
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A  STUDY  OF  EXPLOSIVE  WAVE  PROPAGATION  IN 
GRANULAR  MATERIALS  WITH  MICKOSTRUCTURE 

CHAPTER  1 
INTRODUCTION 

1.1  BACKGROUND 

Traditionally,  in  engineering  practice,  all  stress  analyses  are  conducted 
within  the  framework  of  various  branches  of  continuum  mechanics.  In  doing  30, 
It  is  tacitly  assumed  that  the  microstructural  details  of  the  material  can  be 
neglected.  The  material  is  then  "replaced"  by  an  equivalent  continuum  with 
gross,  or  "overall,"  properties.  The  continuum  approach  has  indeed  been  very 
successful  and  has  led  to  the  development  of  many  useful  theories  of  material 
characterization.  On  the  other  hand,  since  these  theories  disregard  the 
microstructural  details  of  the  materials  under  study,  they  cannot  be  used  to 
determine  how  local  structure  influences  the  gross  behavior  of  the  material. 
This  is  a  real  shortcoming  of  the  continuum  theories,  especially  when  they  are 
applied  for  character ization  of  geological  materials  such  as  sand,  clay  and/or 
rock.  These  types  of  materials  are  commonly  classified  as  materials  with 
microstructure  since,  at  the  micro  level,  the  density,  along  with  other 
Important  variables,  are  not  continuous.  Modeling  of  these  materials  using 
classical  continuum  mechanics  (e.g.,  elasticity,  plasticity,  viscoelasticity, 
etc.)  has  progressed  to  a  point  where  any  fundamentally  new  information  will 
probably  have  to  come  from  a  theory  Incorporating  properties  such  as  grain 
size,  local  porosity,  packing,  etc.  Some  advances  have  been  made  recerttly  in 
developing  analytical  tools  and  models,  which  account  for  some  of  the  struc¬ 
tural  details  of  particulate  materials  such  as  sand.  Two  examples  of  such 
work  are:  (1)  the  so-called  "distributed  body"  concept  advanced  by  Goodman 
and  Cowin  (Reference  31)  and  (2)  "discrete  element"  modeling  pioneered  by 
Cundall  (Reference  20).  The  central  theme  of  the  distributed  body  concept  is 
the  introduction  of  the  "volume  distribution  function"  (a  new  kinematic 
variable)  which  accounts  for  local  porosity  and  its  spatial  gradient.  The 
discrete  element  concept  is  basically  a  numerical  procedure  requiring  large 
computer  simulations  of  grain- to-grain  interaction. 


The  objective  of  this  investigation  was  to  develop  a  theoretical  frame¬ 
work,  and  the  associated  computer  software,  based  on  the  distributed  body 
concept  for  studying  plane  wave  propagation  in  granular  materials  due  to 
airblast  loading.  This  theory  will  allow  specific  relationships  to  be  devel¬ 
oped  between  microstructure  and  wave  propagation  variables.  Wave  propagation 
studies  based  on  the  distributed  body  concept  were  originally  conducted  by 
Nunziato,  Walsh,  et  al .  (References  57-63).  The  present  investigation  will 
extend  their  pioneering  work  to  incorporate  (1)  a  more  realistic  depth- 
dependent  volume  distribution  function  simulating  gravity  effects  in  granular 
soil,  (2)  arbitrary  surface  input  including  finite  times  for  loading  and 
unloading  waves,  and  (3)  probabilistic  considerations  for  treating  non-uniform 
grain  size  and  random  distribution  of  local  porosity. 

1.2  SCOPE 

Chapter  2  contains  a  literature  study  of  previous  work  in  micromechanics 
and  a  summary  of  the  Nunziato,  Walsh,  et  al.,  distributed  body  wave  propaga¬ 
tion  studies.  Extension  of  this  theory  to  include  items  1  through  3  in  the 
above  paragraph  is  documented  in  Chapter  3.  Parametric  results  from  the 
extended  theory  are  presented  in  Chapter  if.  A  summary  and  recommendations  are 
given  in  Chapter  5.  Finally,  a  user's  guide  for  the  microstructural  wave 
propagation  code  MIC1D  is  given  in  Appendix  A. 


CHAPTER  2 

REVIEW  OF  PREVIOUS  WORK 

2.1  GENERAL  BACKGROUND 

Studies  of  geological  materials  with  microstructure  started  many  years 
ago  with  research  on  granular  materials  modeled  as  aggregate  assemblies  of 
discs  or  spheres.  An  excellent  review  article  by  Derestewtcz  (Reference  22) 
presents  both  static  and  dynamic  studtes  prior  to  1958.  Another  more  recent 
review  article  by  Krlzek  (Reference  48)  appeared  In  1971,  and  presented 
basically  the  dynamic  response  of  coheslonless  granular  soils.  Three  recent 
symposia  on  this  subject  (References  16,  45,  and  83)  have  Indicated  renewed 
research  Interest. 

The  concept  of  modeling  granular  media  as  an  array  of  elastic  particles 
(spheres  or  discs)  lead  to  the  Initial  attempts  at  predicting  wave  propagation 
through  such  media.  Early  work  by  Ilda  (References  41  and  42),  Takahashl  and 
Sato  (Reference  79),  Hughes  and  Cross  (Reference  39),  Hughes  and  Kelly 
(Reference  40),  Gassman  (Reference  33)  and  Brandt  (Reference  4)  employed  a 
normal  granular  contact  force  concept.  This  Initial  work  Investigated  the 
propagation  velocity  as  a  function  of  confining  pressure,  particle  size  and 
aggregate  geometry. 

It  was  discovered,  however,  that  the  classical  theory  of  contact  due  only 
to  normal  forces  does  not,  In  general,  accurately  model  real  materials.  With 
this  In  mind,  Duffy  and  Mlndlln  (Reference  26)  proposed  a  theory  for  granular 
media  which  Included  both  normal  and  tangential  contact  forces.  This  theory 
produced  a  nonlinear  and  Inelastic  stress-strain  relation.  Hendron 
(Reference  34)  has  also  done  work  In  this  area. 

More  recent  theories  of  granular  media  have  Included  statistical- 
stochastic  approaches,  e.g.,  Hudson  (Reference  38),  Fletcher  (Reference  29), 

Fu  (Reference  30),  Chambre  (Reference  9),  Varadan,  et  al.,  (Reference  81)  and 
Endley  and  Peyrot  (Reference  28).  Quite  recently,  Mroz  (Reference  53)  and  Kuo 
(Reference  49)  employed  continuum  plasticity  concepts  and  general  contact 
theory  In  an  attempt  to  unify  the  treatment  of  granular  materials  at  both  the 
particulate  and  continuum  levels.  Cundall,  et  al.,  (References  20  and  21) 
proposed  a  numerical  method  called  the  discrete  element  technique  for  granular 
and  rock  assemblies,  and  Brown,  et  al.,  (Reference  5)  have  used  this  approach 
for  rubble  screens.  Morland  (Reference  52)  considered  a  rock/ granular  media 


as  a  regularly  jointed  media  and  used  an  anisotropic  elasticity  approach. 
Particulate  media  has  also  been  studied  by  Hill  and  Harr  (Reference  36)  based 
upon  a  diffusion  equation  derived  from  probabilistic  models,  while  Soo 
(Reference  78)  has  considered  the  dynamic  interactions  of  granules.  Rohanl , 
et  al.t  (References  43,  71,  and  72)  have  been  doing  wave  propagation  research 
in  this  area  for  granular  sands  and  layered  soils  using  continuum  models. 
Endochronic  theories  have  also  been  applied  to  granular  soils,  e.g.,  Read  and 
Valanis  (Reference  70),  Lin  and  Wu  (Reference  50),  and  Bazant,  et  al . , 
(Reference  2).  Studies  have  been  made  of  the  propagation  of  waves  through 
elastic  materials  containing  spherical  inclusions,  e.g.,  Mai  and  Bose  (Refer¬ 
ence  51).  Bleich,  et  al.,  (Reference  3)  employed  an  el as tic- pi  as tic  constitu¬ 
tive  law  to  model  a  specific  geomechanics  boundary  value  problem.  Fluid 
saturated  granular  media  have  been  studied  by  Garg,  et  a.l .  (Reference  32), 
Hsieh  and  Yew  (Reference  37),  Vardoulakis  and  Beskos  (Reference  82)  and 
Zienkiewcz  and  Shiomi  (Reference  84).  Nachllnger  and  Nunztato  (Reference  55) 
used  an  internal  state  variable  approach  to  wave  propagation  problems.  Modern 
mixture  theories  (References  25,  64-66,  68,  and  69)  also  show  some  promise  of 
modeling  porous  and/or  granular  media. 

With  regard  to  experimental  work,  the  method  of  photoelasticity  has  been 
used.  This  particular  method  is  quite  well  suited  for  studying  the  detailed 
load  transfer  between  individual  granules  as  whole  field  data  are  obtained 
during  the  experiment.  Photoelasticity  has  been  used  for  granular  media  by 
Drescher  and  de  Josselin  de  Jong  (Reference  23),  Drescher  (Reference  2*0,  and 
Durelli  and  Wu  (Reference  27).  This  work  was,  however,  only  for  static 
behavior.  The  only  dynamic  analysis  of  granular  media  employing  photoelastic- 
ity  was  performed  by  Rossmanlth  and  Shukla  (Reference  75).  Their  technique 
employed  the  use  of  high  speed  photography  to  record  wave  propagation  through 
an  assembly  of  birefringent  discs. 

Of  the  previous  work,  three  new  constitutive  theories  which  show  special 
promise  in  modeling  granular  and  porous  media  are:  the  so-called  "pore- 
collapse"  models  (References  6-8,  33,  46,  and  47);  the  microstructural  models 
based  upon  "fabric  tensors"  (References  11,  44,  45,  56,  and  67);  and  the 
Goodman-Cowln  distributed  body  approach  (References  1,  12-18,  31,  and  80). 

The  pore-collapse  model,  originally  developed  by  Carrol  and  Holt,  is  based 
upon  the  collapse  of  a  single  pore  within  the  media.  Researchers  at  the 
Sandia  National  Laboratories  have  used  this  approach  with  some  success  to 


model  the  dynamic  response  of  porous  and  granular  media.  This  theory, 
however,  cannot  relate  the  effects  of  neighboring  pores  on  one  another,  and, 
hence,  the  effect  of  pore  distribution  cannot  be  accounted  for.  With  regard 
to  the  Goodman-Cowin  distributed  body  theory,  the  medium  is  assumed  to  be 
distributed  in  space  by  an  independent  kinematical  function  called  the  volume 
distribution  function.  Nunzlato,  Walsh,  et  al . ,  (References  57-63)  have 
applied  this  theory  to  several  wave  propagation  studies  and  found  success  in 
modeling  particular  situations.  Consequently,  this  particular  theory  looks 
quite  fruitful.  The  fabric  tensor  models  proposed  by  Oda,  Nemat-Nasser ,  et 
al . ,  (Reference  67)  also  look  promising:  however,  their  application  to 
specific  boundary  value  problems  appears  to  be  several  years  away.  At 
present,  they  are  looking  at  the  details  of  the  microstructural  fabric,  and 
they  eventually  may  give  insight  as  to  the  nature  of  the  volume  distribution 
function  for  a  Goodman-Cowin  body. 

2.2  DISTRIBUTED  BODY  THEORY 

The  distributed  body  theory  originally  developed  by  Cowin  and  Goodman  was 
constructed  to  allow  a  continuum  theory  to  be  applied  to  materials  with 
noncontlnuous  fields  of  mass  density,  stress,  body  force,  etc.  Thus,  the 
model  could  be  used  to  describe  the  behavior  of  a  wide  variety  of  materials 
having  granular  and/or  porous  structures.  Fundamental  to  the  theory  is  the 
assumption  that,  at  any  point  in  the  material,  the  overall  mass  density  p 
may  be  written  as 

p  -  vY  (2.1) 

where  y  is  the  density  of  the  granules  (or  matrix  material)  and  v  *  v  (X,t) 
is  referred  to  as  the  volume  distribution  function.  This  function  describes 
the  way  the  medium  is  distributed  in  space  allowing  for  voids  or  other 
particular  granular  structures.  Thus,  the  theory  uncouples  the  mass  density 
of  the  granules  from  the  mass  density  of  the  entire  material,  and  allows 
compressibility  due  to  both  granule  compressibility  and  void  compaction.  In 
general ,  0  <  v  <  1  ,  and  v  is  related  to  the  porosity  n  and  void  ratio 


Within  a  one- dimensional  framework,  the  classical  balance  law  of  conser¬ 
vation  of  linear  momentum  reads 

pox=3X+pob  (2*3) 

where  T  is  the  stress,  b  is  the  body  force,  x  is  the  particle  position, 
x  is  the  particle  acceleration,  X  is  the  reference  position  coordinate,  and 
subscript  o  denotes  values  in  reference  state.  In  addition  to  thl3 
classical  balance  law,  the  distributed  body  theory  also  requires  an  indepen¬ 
dent  balance  equation  governing  the  volume  distribution.  In  one  dimension, 
this  second  equation  governing  void  change  is  given  by 

pokv“H+pog  f2-** 

where  k  is  called  the  equilibrated  inertia,  h  the  equilibrated  stress  and 
g  the  Intrinsic  body  force.  Physical  Interpretation  of  the  mlcrostruetural 
variables  k  ,  h  and  g  is  somewhat  difficult  to  make.  In  general,  these 
variables  are  related  to  the  local  contact  mechanics  at  the  granular  level  and 
can  be  related  to  particular  self-equilibrated  singular  stress  states  from 
classical  elasticity  (e.g.,  double  force  systems,  centers  of  dilatation).  It 
has  been  proposed  (Reference  63)  that  k  is  related  to  the  void  mean  surface 
area  and  to  the  number  of  voids  present,  h  is  a  result  of  the  interaction 
forces  between  nieghboring  voids  and  will  vanish  when  the  voids  are  suffi¬ 
ciently  separated,  and  g  is  related  to  the  coupling  between  the  total 
deformation  of  the  medium  and  the  changes  in  void  volume. 

For  granular  geological  materials,  we  assume  that  the  media  is  composed 
of  compressible  granules  at  relatively  high  confining  pressures  so  as  to 
prevent  material  flow.  For  this  case,  an  appropriate  constitutive  formulation 
would  read 

T  -  T  <v0,  v,  e)  (2.5) 

and,  hence,  the  stress  depends  upon  the  reference  and  current  volume  distri¬ 
butions,  the  gradient  of  the  volume  distribution,  and  the  strain  e  .  An 
explicit  form  of  Equation  2.5  which  has  been  proposed  (References  31  and  61) 
uses  an  even  quadratic  form  in  the  gradient  of  v  , 


(2.6) 


T  -  v  £A(v0,  v,  e)  ♦  j  a(vQ,  v*  3 

where  A  and  a  are  two  material  functions  of  the  indicated  variables. 
First  and  second  order  moduli  defined  by 
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will  be  needed  for  subsequent  wave  analyses.  Normally  E  >  0  ,  but  the  second 
order  modulus  E  may  be  positive  or  negative. 

2.3  WAVE  PROPAGATION  WITHIN  A  DISTRIBUTED  BODY 

As  previously  mentioned,  the  wave  propagation  theories  set  forth  by 
Nunzlato,  Walsh  and  coworkers  for  Goodman-Cowtn  distributed  bodies  appear  to 
have  excellent  promise  for  application  to  granular  geological  materials.  This 
section  will  briefly  review  some  basic  details  and  previous  results  about 
these  theories. 

The  basic  premise  of  this  particular  wave  theory  lies  In  modeling  the 
wave  as  a  propagating  singular  surface  across  which  there  exists  a  Jump 
discontinuity  In  a  particular  variable.  Commonly  dynamic  loadings  will 
produce  second-order  acceleration  waves,  having  a  jump  discontinuity  In  the 
particle  acceleration  at  the  wave  front.  In  some  cases,  however,  the  loading 
could  produce  a  fir3t-order  shock  wave,  having  a  jump  In  the  particle  velocity 
at  the  wave  front.  Most  modeling  In  these  materials  has  been  done  for  accel¬ 
eration  waves,  and  this  case  will  now  be  described. 

As  mentioned,  a  wave  is  modeled  as  a  propagating  singular  surface  of  zero 
thickness  with  speed  U  ,  see  Figure  2.1.  The  jump  of  a  quantity  <|>  across 
this  surface  Is  defined  by 


[*]  -  ♦'  -  *+  (2.8) 

where  $+  and  are  the  limiting  values  of  $  Immediately  ahead  of  and 

behind  the  wave.  An  acceleration  wave  Is  therefore  defined  as  a  wave  across 
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which  the  particle  velocity,  strain,  and  volume  distribution  are  continuous 
but  their  spatial  and  temporal  derivatives  are  not.  Thus,  this  type  of  motion 
carries  propagating  discontinuities  in  the  particle  acceleration  and  various 
other  gradients  of  the  strain  and  volume  distribution.  The  jump  in  the  parti¬ 
cle  acceleration  [x]  is  called  the  wave  amplitude,  and  will  be  denoted  by 
a(t)  .  Note  that  for  compressive  waves,  a(t)  >0  ,  while  for  expansive 
waves,  a(t)  <  0  . 

Following  singular  surface  analysis  procedures  which  have  now  become 
somewhat  standardized  (Reference  10),  Nunziato,  et  al . ,  developed  the  follow¬ 
ing  expressions  for  two  different  types  of  waves 
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with  subscripts  e  ,  v  ,  and  X  meaning  partial  differentiation  with  respect 
to  the  Indicated  variable,  and  (•)  meaning  immediately  ahead  of  the  wave. 
The  speed  UF  denotes  the  "fast"  wave  speed  which  is  associated  predominantly 
with  the  elasticity  of  the  granules.  The  quantity  Ug  is  the  "slow"  wave 
speed  and  is  connected  to  the  compressibility  of  the  material  due  to  consoli¬ 
dation. 

The  wave  amplitude  a  ,  which  is  equal  to  the  jump  discontinuity  in  the 
acceleration  across  the  wave  front,  has  also  been  studied.  Nunziato,  et  al . , 
have  found  that  the  amplitude  for  one-dimensional  wave  propagation  satisfies 
the  following  nonlinear  Bernoulli  equation 


-  K(X)a2  -  p(X)a 


(2.11) 


where  y(X)  and  k(X)  are  material  coefficients  given  in  general  by  rather 
lengthy  expressions.  The  coefficient  p(X)  is  related  to  dispersive  effects, 


*  .  ■■ 


while  <(X)  reflects  both  the  elastic  response  of  the  granules  and  dispersive 
effects.  Depending  on  the  nature  of  <  and  y  ,  the  theory  can  predict 
growth  or  decay  of  wave  amplitude. 

Nunzlato,  et  al.,  (Reference  61)  presented  an  application  of  these 
general  theoretical  results  to  a  specific  granular  medium,  PBX-9404  (an  explo¬ 
sive  powder).  They  chose  a  volume  distribution  to  be  periodic  in  nature. 


l.e., 


V0<x>  -  va  ♦  (1  -  va>  COS^ 


(2.12) 


where  v  is  a  material  constant,  and  l  is  a  characteristic  length  presuma- 
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bly  related  to  the  grain  size.  For  PBX-9404,  v  -  0.984  and 

3 


t  »  1.5  mm  were  chosen.  The  specific  constitutive  form  for  this  application 
was  selected  as  that  given  in  Equation  2.6. 

Using  the  previous  specific  forms  and  assuming  that  the  wave  starts  at 
X  -  0  in  a  granule,  the  fast  wave  speed  becomes 
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where  Ug  is  the  wave  speed  in  a  granule  and 
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In  attempting  to  compare  with  experimental  data,  Nunziato,  et  al., 
(Reference  61)  point  out  that  what  is  actually  measured  is  the  transit  time  of 
the  wave  t  .  This  quantity  is  a  function  of  the  propagation  distance  X  , 
and  l 3  related  to  the  average  wave  velocity  U  by  the  expression 


(2.15) 


Applying  this  to  the  fast  wave,  one  can  write 
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which  can  be  expressed  In  terms  of  an  elliptic  Integral  of  the  first  kind  F  , 
l  .e., 

t(X)  "  2^0“  F<2lfX/*>M)  (2-17) 

g 

Using  the  boundary  condition  t(4)  =  i  ,  with  known  values  of  U.  ,  t  ,  and 

S  ®  g 

l  yields  an  equation  whose  root  gives  the  value  of  M  .  For  PBX-9404,  with 

U  «  3.71  km/s  and  t  -  0.49311  ys  ,  Equation  2.17  yields  M  -  0.7525  . 

8  g 

The  amplitude  behavior  Is  governed  by  Equation  2.11  and,  for  this 
specific  case,  the  coefficients  <  and  y  become 


Again,  for  the  specific  material  PBX-9404,  Nunzlato,  et  al.,  using  low- 

amplitude  shock  wave  experiments,  found  that  (A  )  =  -58  GPa  and 

s  2  ee  o 

(oee)0  -  7.31  x  10  GPa  -  mm  . 

Further  developments  of  this  theory,  along  with  a  general  purpose  compu¬ 
ter  program,  have  been  developed  by  Sadd  (Reference  76),  to  evaluate  the 
average  wave  speed  and  amplitude  behavior.  Typical  results  for  the  specific 
material  values  Tor  PBX-9404  are  shown  In  Figures  2.2  and  2.3.  Figure  2.2 
Illustrates  the  behavior  of  the  average  wave  speed  with  propagation  distance. 
It  Is  evident  that  the  mlcrostructural  effects  predominate  at  Initial  dis¬ 
tances  producing  a  large  variation  In  wave  speed.  Gradually,  as  the  wave 
moves  further  Into  the  medium,  the  speed  has  less  variation  and  approaches  a 

constant  value.  The  amplitude  behavior  with  propagation  distance  13  shown  In 

2 

Figure  2.3  for  the  case  of  an  Initial  amplitude  of  2.7  Gm/s  .  For  this 
case,  the  amplitude  decays  with  a  superimposed  periodic  oscillation. 

Nunzlato,  et  al . ,  (Reference  61)  also  presented  actual  speed  and  amplitude 
experimental  data  for  this  material  and  found  fairly  good  agreement  with  the 
theory . 


SINGULAR  SURFACE  WAVE 


Propagating  Discontinuity:  [«()]*  <J>  —  4»"l_ 

Wave  Amplitude:  a  =  [x] 


Figure  2.1.  Schematic  of  a  propagating  singular  surface. 
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CHAPTER  3 


i 

!  DEVELOPMENT  OF  WAVE  PROPAGATION  THEORY 


3.1  GENERAL 

The  major  purpose  of  the  research,*  herein  reported  was  to  develop  a  one- 
dtmenstonal  wave  propagation  theory  and  associated  computer  code  incorporating 
the  distributed  body  theory  to  account  for  material  microstructure.  The 
developed  theory  is  general  in  that  it  can  handle  a  variety  of  volume  distri¬ 
bution  functions,  and,  thus,  it  can  model  several  types  of  microstructure.  In 
addition,  the  theory  has  been  extended  to  include  the  wave  motion  description 
of  more  than  simply  one  singular  surface  as  discussed  previously  In 
Section  2.3.  In  this  regard,  a  wave  profile  constructed  of  several  singular 
surface  waves  has  been  analyzed.  An  uncoupled  theory  In  which  each  wave 
propagates  independently  has  been  developed.  In  addition,  a  coupled  theory 
incorporating  microstructural  changes  with  the  passage  of  each  wave  in  a 
profile  has  also  been  constructed.  The.finlal  step  of  the  theoretical  work 
was  to  construct  a  probabilistic  analysis  based  upon  the  developed  code  using 
a  moment-generating  procedure  due  to  Rosenblueth  (Reference  73).  The  probabi¬ 
listic  analysis  allows  for  the  treatment  of  grain  size  and  local  porosity  as 
random  variables.  The  following  sections  discuss  in  detail  each  of .these 
developments. 

3.2  VOLUME  DISTRIBUTION  FUNCTIONS 

In  order  to  apply  the  distributed  body  theory  and  develop  a  wave  propaga¬ 
tion  analysis,  it  is  necessary  to  have  explicit  constitutive  forms,  see  for 
example  Equation  2.6,  and  the  initial  volume  distribution  \>Q(X)  must  also  be 
specified.  Any  proposed  volume  distribution  function  should  reflect  the 
density  variations  and  other  microstructural  features  within  the  material.  It 
i3  difficult  to  construct  such  a  function  which  characteri zes  these  variations 
precisely  and  yet  has  the  smoothness  requirements  to  be  compatible  with  the 
theory.  We  will  follow  the  approach  that  vQ(X)  should  be  a  continuous  func¬ 
tion  in  order  to  perform  certain  required  differentiations  and  integrations 
and  that  it  yield  the  correct  average  density. 

As  discussed  previously  in  Section  2.3.  Nunzlato,  et  al . ,  (Reference  61) 
in  constructing  a  wave  propagation  analysis,  developed  a  specific  volume 
distribution  function.  Their  work  was  for  a  granular  material,  PBX-9404,  an 


explosive  powder*/ binder  system.  They  proposed  a  periodic  structure  of  the 
form 

voOO  =  va  ♦  (1  -  va)  cos  ^  (3.1) 

where  v  and  l  are  material  constants.  This  function  Is  plotted  In 

el 

Figure  3.1  versus  the  distance  coordinate. 

The  quantity  v  would  be  given  by  the  overall  density  of  the  material 
divided  by  the  granule  density  and  Is  thus  related  to  the  average  value  of  the 
volume  distribution.  The  second  material  constant  *,  Is  referred  to  as  a 
characteristic  length  associated  with  this  periodic  structure.  Clearly  l 
specifies  the  length  of  the  repeating  units  of  the  microstructure.  For 
granular  materials,  l  would  be  related,  but  not  necessarily  equal,  to  the 
average  grain  size. 

In  regard  to  this  characteristic  length,  the  work  of  Shahlnpoor  (Refer¬ 
ence  77)  Is  appropriate  to  consider.  Shahlnpoor  did  experiments  of  randomly 
packed  spherical  granules  on  a  flat  surface.  His  work  demonstrated  the 
concept  of  distinct  packing  geometries  referred  to  as  Voronol  cells,  see 
Figure  3.2.  It  Is  evident  that  for  some  packing  geometries,  If  a  periodic 
structure  Is  assumed,  the  characteristic  length  Jl  ,  being  equal  to  the 
Voronol  cell  size,  could  be  several  grain  diameters. 

Since  the  mechanical  response  of  most  geological  materials  like  sand  or 
gravel  Is  affected  by  in  situ  conditions  such  as  overburden,  the  microstruc¬ 
ture  will  be  nonhomogeneous ,  l.e.,  be  depth  dependent.  With  this  in  mind,  a 
new  volume  distribution  function  was  developed  which  can  predict  such  a  struc¬ 
ture.  One  particular  form  uses  an  exponential  factor  and  may  be  written  as 

vQ(X)  =  1  -  (1  -  vb)e"BYX  (3.2) 

where  vb  ,  B  ,  and  y  are  material  constants.  A  plot  of  this  distribution 
function  is  shown  in  Figure  3.3.  Clearly  for  this  case,  the  material  becomes 
more  dense  with  depth  X  into  the  medium.  The  constant  is  the  volume 

distribution  at  the  free  surface  X  *  0  ,  y  is  the  average  density  of  the 
material,  and  the  constant  B  determines  the  rate  of  consolidation  with 
depth.  It  3hould  be  pointed  out  that  this  exponential  form  does  not  contain 


any  periodic  structure;  hence,  it  should  produce  monotonic  results  for  the 
wave  propagation  characteristics . 

Another  volume  distribution  function  which  was  used  involves  the  combina¬ 
tion  of  the  periodic  form  given  by  Equation  3.1  and  the  exponential  form  from 
Equation  3.2.  The  combined  form  involves  simply  the  product  of  these  two 
relations,  i.e., 

voOO  =  [va  ♦  (1  -  va)  cos  ^p][1  -  (1  -  vb)  e'BYX]  (3.3) 

and  again  v  ,  \>.  ,  g,  ,  B  and  Y  are  matertal  constants.  It  is  evident 

a  D 

that  this  form  (shown  in  Figure  3.H)  will  thus  produce  a  combined  periodic- 
exponential  depth  dependent  microstructure. 

During  the  course  of  this  investigation,  other  forms  of  the  volume 
distribution  function  were  developed  Including  algebraic  and  additive 
per  iodic- exponential  forms.  However,  the  three  forms  given  by  Equations  3.1- 
3.3  appear  to  provide  a  broad  enough  microstructure  model  for  the  objectives 
of  this  research.  Consequently,  only  these  three  form 3  will  be  included  in 
the  remaining  sections  of  this  report. 

3.3  WAVE  PROPAGATION  ANALYSIS 

The  wave  propagation  analysis  and  the  development  of  an  associated 

computer  code  was  done  based  upon  the  previous  fundamental  work  of  Goodman  and 

Cowin  (Reference  31)  and  Nunziato,  et  al .  (Reference  61  ) .  The  constitutive 

form  given  by  Equation  2.6  was  also  used  in  this  work.  Equation  2.6  was  used 

by  Nunziato,  et  al.,  but  was  originally  proposed  by  Goodman  and  Cowin  in  1972. 

3v 

The  constitutive  dependence  on  the  gradient  of  the  volume  distribution  is 

3  v 

significant  and  allows  an  equilibrium  stress  to  depend  on  -rrr  .  Since  Equa- 

7k 

tion  2*6  Involves  the  square  of  ^  ,  It  will  be  an  isotropic  form  in  that 

o  A 

variable  (required  by  material  frame  indifference)  and,  hence,  the  stress 
response  will  be  independent  of  the  sign  of  the  gradient.  Also,  the  presence 
of  the  gradient  term  allows  the  theory  to  predict  a  generalized  Mohr-Coulomb 
failure  criterion  (Reference  31). 

Obviously,  the  two  material  functions  A  and  a  defined  in  Equation  2.6 
will  specify  the  response  of  the  medium  to  deformation.  Equation  2.6  indi¬ 
cates  that  the  material  function  a  ,  specifies  the  effect  of  the  gradient  of 
the  volume  distribution.  If  a  is  small,  then  the  stress  will  not  be 
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significantly  influenced  by 
strain  state), 


3v 

3X  * 


For  a  stress-free  reference  state  (zero 


A(v  .  vQ  ,  0)  =  0 

o(v0,  vQ  ,  0)  =  0  (3.*0 

Considering  the  stress-strain  behavior  which  could  come  from  Equa¬ 
tion  2.6,  Figure  3.5  illustrates  some  typical  curves  for  various  volume 
distribution  functions.  It  should  be  pointed  out  that  the  shape  of  this  curve 
could  vary  considerably  for  various  types  of  geological  materials  and  is  a 
function  of  rate  of  loading.  This  figure  demonstrates  the  stress-strain 

behavior  of  the  granular  assembly  medium  accounting  for  the  particular 

3  v 

reference  values  of  v(X)  and  ,  i.e., 

T  -  T(X,e)  =  T(vo,  vQ,  (|^)Q,  e)  (3-5) 

From  such  typical  behavior  as  shown  in  Figure  3-5,  it  is  apparent  that 
the  two  moduli  E  and  E  given  by  Equation  2.7  would  satisfy  the  relations 


e-£>o 


3e 


(3-6) 


In  particular,  from  the  theory,  it  can  be  shown  that  the  fast  wave  speed 
given  by  Equation  2.9,  can  also  be  written  as 


UF  -Vvv7  0.7) 

where  Yq  is  the  density  of  the  granules,  and  E  is  to  be  evaluated  in  the 
reference  3tate.  Consequently,  for  real  wave  speeds  E  >  0  ,  and  from  Equa¬ 
tion  2.7,  this  means  that 


where  subscripts  on  A  and  a  mean  partial  differentiation.  Equation  3.9 
then  gives  a  condition  that  the  material  functions  A  and  a  must  satisfy 
for  a  given  volume  distribution. 

With  regard  to  the  amplitude  behavior,  It  was  shown  that  the  amplitude 
obeyed  the  differential  equation  given  by  Equation  2.11.  By  combining  Equa¬ 
tions  2.7  with  2.18,  it  becomes  apparent  that  the  coefficient  k  Is  related 
to  E  by 


K 


2Vf 


(3.10) 


and,  hence,  the  curvature  of  the  stress-strain  curve  will  affect  the  amplitude 
behavior. 

In  order  to  further  elaborate  on  the  constitutive  relationship  of  this 
theory,  consider  the  following  special  cases: 

(1 )  Uniform  Volume  Distribution: 

For  this  case,  v  =  v  =  constant,  and  so  Equation  2.6  yields 

T  =  T(e)  »  vA(v,  v,  e)  (3.11) 


which  can  be  Interpreted  as  a  nonlinear  elastic/ plastic  material.  Conse¬ 
quently,  the  material  parameter  A  is  associated  with  constitutive  behavior 
of  the  microstructural  media  based  upon  the  local  volume  distribution  but 
neglecting  distribution  gradients.  Wave  propagation  studies  for  this  case 
reduce  to  the  classical  one-dimensional  plastic  wave  motion  analyses  (see 
Cristescu,  Reference  19). 

(2)  Homogeneous  Elastic  Case: 

For  the  reduction  to  linear  elasticity,  the  volume  distribution  is  taken 
v  *>  1  .  Hence,  from  Equation  3.11 


to  be  unity,  l.e., 


and  for  the  linear  elastic  case 


A(1,  1,  e)  -  Ee  (3.13) 

where  E  Is  the  elastic  modulus  and  Is  Identical  to  the  first  order  modulus 
previously  defined  In  Equation  2.7.  Wave  motion  analysis  for  this  case  yields 
the  well  known  results 

u  =  VI7p 

a  =  aQ  =  constant  (3. 11*) 

Based  upon  this  wave  motion  analysis,  a  computer  code*  was  developed  to 
handle  any  of  three  volume  distribution  functions  given  by  Equations  3. 1-3-3- 
The  constitutive  form  Incorporates  Equation  2.6,  with  specific  values  for  the 
two  material  functions  A  and  a  to  be  Input  by  the  user.  The  code  uses 
general  techniques  of  numerical  Integration  using  four-point  Gauss  quadrature 
to  calculate  the  necessary  Integrals  for  computation  of  the  average  wave 
speed,  see  Equations  2.15  and  2.16.  In  addition,  a  fourth  order  Runge-Kutta 
scheme  Is  used  to  solve  the  nonlinear  amplitude  Equation  2.11.  Thus,  the 
basic  features  of  the  code  were  to  calculate  the  wave  speed  and  amplitude 
(particle  acceleration)  at  various  positions  and  times. 

Typical  results  of  the  code  are  shown  In  Figures  3.6-3. 14.  The  first  set 
of  figures  (Figures  3-6~3-8) ,  Illustrates  results  using  the  periodic  volume 
distribution  function  given  by  Equation  3.1.  Recall  this  distribution  func¬ 
tion  was  shown  in  Figure  3.1.  The  specific  material  parameters  for  these 

_  %  2  <t 

results  are  va  -  0.85  ,  l  »  0.1  in  ,  y  *  2.4  x  10  lb-sec  /in  f 
a  o 

9  5  2 

a  »  -450  lb  ,  a  -  1  x  1 0  lb  ,  A  -  3  x  1 0  lb/in  , 
e  ee  e 

• 

A  -  -1  X  10  lb/ln*  .  This  material  will  be  referred  to  as  Material  PI. 
ee 

The  actual  wave  speed  behavior  shown  In  Figure  3.6  varies  periodically,  as 
given  by  Equation  2.13.  However,  the  average  wave  speed  will  oscillate 
Initially  and  then  approach  a  constant  value,  as  shown  in  Figure  3.7.  The 
behavior  of  the  amplitude  ratio  (normalized  by  ao  )  Is  shown  in  Figure  3.8 


*  The  computer  code  Is  referred  to  as  MIC1D. 


and  Illustrates  the  expected  effect  of  the  Initial  amplitude  on  the  attenua¬ 
tion  rate;  l.e.,  the  higher  aQ  ,  the  larger  the  attenuation. 

Figures  3*9-3.11  show  the  corresponding  results  for  the  exponential 
volume  distribution  specified  In  Equation  3.2  and  shown  In  Figure  3.3. 

Material  parameters  for  this  case  are  -  0.65  , 

_il  *  H  C  • 

If  -  2.4  x  10  lb-sec  /In  ,  a  -  -8  x  10  lb  ,  o  -  4.8  x  10  lb  , 
o  e  ee 

3  2  2  2 

A  -  3  x  10  lb/in  ,  A  -  -1.0  lb/ln  ,  B  -  10  In  /lb  , 
e  ee  ’ 

-2  * 

Y  -  7.2  x  10  lb/ln  .  This  material  will  be  referred  to  as  Material  El. 

The  wave  speed  (shown  in  Figures  3.9  and  3. 10)  is  now  a  monotonlcally 

Increasing  function  with  depth  X  since  the  porosity  Is  decreasing  in  that 
direction.  Furthermore,  the  amplitude  behavior  shown  in  Figure  3. 11  Illus¬ 
trates  a  much  less  pronounced  attenuation  rate  when  compared  with  the  periodic 
volume  distribution  results  In  Figure  3.8.  The  reason  for  this  behavior  is 
the  fact  that,  for  the  exponential  distribution  function,  the  material 
response  rapidly  approaches  with  depth  that  of  an  elastic  material. 

Finally,  results  of  using  the  combined  periodic-exponential  volume 
distribution  function  given  by  Equation  3.3  are  shown  in  Figures  3.12-3.14. 
Again,  this  particular  distribution  function  was  shown  previously  in 
Figure  3.4.  Model  parameters  for  this  case  are  v  *  0.992  ,  v.  -  0.65  , 

3  D 

£,  -  0.06  In  ,  YQ  -  2.4  x  10~4  lb-secVln”  ,  Y  -  7.2  x  10~2  lb/ln’  , 

B  -  30  In’/lb  ,  a  -  -3  x  10*  lb  ,  a  -  5  x  10*  lb  ,  A  »  3  x  10*  lb/ln*  , 

e  ee  e 

2 

A  -  -750  lb/in  .  This  material  will  be  referred  to  as  Material  PEI. 
ee 

Results  for  wave  speed  and  amplitude  attenuation  indicate  combined  features  of 
each  of  the  two  previous  distribution  functions. 

Additional  features  to  calculate  particle  velocity  and  displacement, 
stress,  wave  profile  behavior,  and  probabilistic  effects  were  also  added  to 
the  basic  code.  These  developments  are  discussed  In  the  next  two  sections. 

3.4  WAVE  PROFILE  ANALYSIS 

This  section  describes  the  efforts  to  extend  the  basic  theory  to  predict 
wave  profile  behavior  where  the  wave  would  have  a  definite  rise  time.  This 
situation  requires  that  consideration  be  given  to  a  train  of  waves  moving 
together.  The  previous  modeling  of  treating  a  wave  as  a  singular  surface  of 
zero  thickness  and  duration  must  be  modified.  Figure  3. 15  illustrates  the 


procedure  of  constructing  a  profile  from  a  series  of  Impulsive  singular  waves. 
The  central  complication  In  this  procedure  Is  the  fact  that  wave  analyses  are 
required  for  the  cases  of  waves  traveling  behind  one  another.  What  this  means 
Is  that,  with  the  exception  of  the  leading  wave,  all  waves  will  be  moving  Into 
media  which  are  undergoing  non- stationary  deformation.  If  the  wave  Is  moving 
Into  a  region  which  Is  not  at  rest  In  Its  reference  configuration,  then  the 
analysis  for  the  wave  speed  and  the  amplitude  attenuation  will  be  greatly 
complicated. 

A  simple  example  of  this  complication  may  be  seen  from  the  wave  speed 
relations  given  In  Equation  2.9.  The  velocity  of  propagation  was  given  by  an 
equation  containing  terms 
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e  -  (J-)  — r~ 
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(3.15) 


with  T 


and  T^  being  moduli  evaluated  Immediately  in  front  of  the 


given  wave  and  where  T 


,  and  T 


Clearly,  the 


state  of  the  material  ahead  of  the  wave  as  specified  by  the  terms  T£  ,  h^ 
T+  will  have  a  complicating  effect  on  the  calculation  of  the  wave  speed. 


Note  that  the  quantity  g  will  vanish  if  the  wave  Is  moving  into  a  region 
which  is  stress  free.  The  state  of  affairs  is  considerably  worse  for  the  case 
of  the  wave  amplitude  analysis  where  the  coefficients  of  Equation  2.11  become 
quite  long  and  complicated  functions  of  the  deformation  state  In  front  of  the 
wave. 

It  was  decided  that,  In  light  of  the  time  restrictions  of  the  current 
investigation,  the  analysis  of  constructing  a  profile  from  a  group  of  travel¬ 
ing  waves  be  made  under  the  simplifying  assumption  that  the  propagatlonal 

characteristics  of  each  wave  depend  solely  upon  the  volume  distribution  at  the 

« 

wave  front.  This  volume  distribution,  In  turn,  depends  upon  the  current 
stress  state  at  the  wave  front.  What  this  means  is  that  waves  traveling 
behind  the  leading  wave  will  feel  a  different  material  caused  by  the  change  of 
stress  due  to  all  previous  waves. 


In  order  to  employ  this  modeling  concept,  a  method  to  compute  the  stress 
associated  with  a  singular  surface  wave  must  be  developed.  Recall  that  the 
amplitude  of  the  wave  was  originally  defined  as  the  jump  in  the  media  acceler¬ 
ation,  l.e., 


a  -  [x]  »  x  -  x 


(3.16) 


where  x  and  x  are  the  limiting  values  of  the  acceleration  just  ahead  of 
and  behind  the  singular  surface  wave.  The  equation  of  motion  was  given  by 


V  ■  SX  *  pob 


(3.17) 


where  T  is  the  stress  and  b  is  the  body  force  which  is  continuous  every¬ 


where  . 


Using  the  basic  definition  for  Equation  3-16,  we  evaluate  the  jump  of 
Equation  3-17  across  a  typical  wave 


P0[x]  ‘  +  <>oCb] 


which,  if  the  body  force  is  continuous,  can  be  written  as 


v  a  /3T.-  ,3Tx  + 

Voa  “  3X^  (3XJ 


and,  thus,  we  can  write  an  expression  for  the  stress  gradient  behind  a  given 
wave  in  terms  of  the  gradient  in  front  of  the  wave  as 


<I>"  -  Voa  ♦  (i}+ 


(3.18) 


Next,  by  using  a  simple  differencing  scheme 


T  -  T 
n+1  n  , 3T  N- 

(AX)n  "  dX^n 


(3.19) 


and  so, 
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n+1  9X  n  n  n 


(3.20) 
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Finally,  combining  Equations  3*18  and  3.20  gives 


T  i  =■  T  +  Cv  Y  a  +  (||)I3(AX) 
n+ 1  n  oon  9Xn  n 


(3.21) 


Hence,  if  we  know  the  stress  at  one  wave  Tn  ,  we  can  compute  the  new  value 


Tn+1  at  the  next  wave. 


As  an  example  to  Implement  this  theory,  consider  the  M-wave  profile  as 
shown  in  Figure  3.16.  The  stress  values  for  this  case  using  Equation  3.21 


follow  to  be 


T  =  0 
o 


T'  * 


^0TQa i ( AX  j ) 


T2  *  [ v0^0a i  ♦  (||)t](AX2) 


Ta  +  l-v0"^0a2  +  (gjj)  1  3  (AX  2  ) 


-  v0Y0[aj (AX , )  +  (a,  +  a2)(AX2)] 


Ti.  =*  T,  +  tvoYQa,  +  (■^)33(AXs) 


v^CajAXj)  +  (a,  +  a2)  (AX2)  +  (ax  +  a2  +  a,)(AX,)] 


(3.22) 


For  the  general  case  with  n  >  1  ,  the  stress  13  given  by 


Tn  -  v0YQCa, (AX, )  +  (a,  +  a2)(AX2)  +  (a,  ♦  a2  ♦  a,)(AX,) 


,+  (a,  +  a2  +  as+  ...  +  an_1 )(AXn_i )3 


(3.23) 


Now,  since  it  is  expected  that  tf\e  stress  will  affect  the  microstructure, 
we  postulate  that  there  must  be  some  relationship  between  the  average  volume 
distribution  function  ^  and  the  stress  T  ,  l.e.. 


v  -  v(T) 


(3. 24) 


With  little  or  no  stress  v  -  vQ  and  as  the  stress  Increases  (compression 
positive),  one  would  assume  that  v  1  .  Figure  3.17  Illustrates  such 
behavior  for  a  typical  sandy  soli.  Based  upon  these  Ideas,  the  quantity 
v  In  the  periodic  distribution  forms  was  modified  as  a  function  of  the  stress 

di 

by  the  relation 

v  -  1  -  0  "  v  )  e-MT  (3.25) 

a  ao 

where  v  Is  Its  reference  value  and  M  is  a  material  constant.  This 
ao 

simplified  approach  is  essentially  varying  porosity  with  stress  to  predict 
wave  coupling  effects.  This  should  be  regarded  as  an  approximate  technique 
since  Equation  3.25  might  not  be  strictly  compatible  with  the  basic  constitu¬ 
tive  form  (Equation  2.6). 

Equation  3.25  was  then  placed  Into  the  code  to  provide  an  approximate 

means  to  calculate  the  wave  propagatlonal  characteristics  of  waves  traveling 

behind  each  other.  Of  course,  an  uncoupled  theory  would  be  generated  by 

specifying  M  -»  0  which  gives  v  -  v_  and,  thus,  all  waves  will  travel 

a  aQ 

independent  of  one  another.  By  tnputing  a  number  of  waves  of  various  Initial 
amplitudes  with  equal  Initial  time  spaclngs  Ato  ,  a  wave  profile  (accelera¬ 
tion  versus  time)  for  various  depths  can  be  constructed.  With  Equation  3.23, 
the  stress  profile  can  also  be  constructed.  Finally,  from  the  acceleration 
profiles,  the  velocity  and  displacement  profiles  were  calculated.  These 
profile  constructions  are  contained  within  the  code.  Figures  3*18-3.21  Illus¬ 
trate  some  typical  profiles  for  the  periodic  volume  distribution  case  using 
model  parameters  of  Material  PI.  These  results  are  for  the  uncoupled  case 
M  -  0  . 

3.5  PROBABILISTIC  CONSIDERATIONS 

The  purpose  of  a  probabilistic  analysis  Is  to  develop  a  method  by  which 
the  variability  or  uncertainties  In  the  independent  (Input)  parameters  In  a 
particular  problem  can  be  evaluated  or  estimated  in  terms  of  their  effects  on 
the  dispersion  of  the  dependent  (output)  variables.  An  extremely  useful 
procedure  for  determining  the  moments  of  a  dependent  variable  in  terms  of 
functions  of  the  moments  of  Its  Independent  variables  was  developed  by  Rosen- 
blueth  (Reference  73).  The  Rosenblueth  procedure  is  quite  versatile  and  Is 
not  bound  by  the  restrictions  often  Imposed  on  other  moment  generating 


procedures,  such  as  the  method  of  partial  derivatives.  If  Y  Is  functionally 
related  to  two  randan  variables  X!  and  X2  ,  l.e., 

Y  -  Y  (X,  ,  X2)  (3.26) 

and,  If  X!  and  X2  are  uncorrelated  and  their  probability  distribution  func¬ 
tions  are  symmetrical,  then,  according  to  the  Rosenblueth  procedure,  the 
expected  value  of  Y  ,  E(Y)  ,  and  variance  of  Y  ,  V(Y)  ,  can  be  estimated 
from  the  following  expressions 

E(Y)  =  (1/4)  (Y++  +  Y+~  +  Y_+  +  Y  )  (3-27) 

V(Y)  -  E(Y*)  -  [E( Y ) ]*  (3.28) 

where 

E(Y2)  -  (1/4)  [(Y++)2  ♦  (Y+_)2  +  (Y~V  +  (Y-")2]  (3.29) 

and 

Y**  »  Y (X,  ±  oY  ,  X2  ±  aY  )  (3-30) 

A  i  A2 

In  Equation  3*30,  X,  and  X2  are  the  expected  values  of  the  random 
variables  X,  and  X2  ,  respectively.  Slmllary,  oY  and  <j  are  the 
standard  deviations  of  the  random  variables  X,  and  X2  .  It  should  be  noted 
from  Equations  3.27-3.30  that  the  expected  value  and  the  variance  of  Y  can 
be  calculated  from  four  (22)  "point  estimates"  of  the  function  Y  ,  a3  stipu¬ 
lated  by  Equation  3.30.  Each  of  these  point  estimates  can  be  viewed  as  a 
"deterministic  calculation"  using  the  dependent  random  variable  Y  .  The 
above  system  of  equations  can  readily  be  generalized  to  n  random  variables 
requiring  2n  point  estimates. 

A  major  advantage  of  the  Rosenblueth  procedure  over  the  Monte  Carlo 
method  can  now  be  realized  when  comparing  the  minuter  of  deterministic  calcula¬ 
tions  required  to  determine  the  moments  of  the  random  variable  Y  .  For 

example,  In  the  case  of  three  random  variables,  the  Rosenblueth  procedure 

$ 

requires  2  -  8  calculations.  The  Monte  Carlo  method,  on  the  other  hand, 

may  require  several  hundred  calculations.  It  should  also  be  pointed  out  that 


the  Rosenblueth  procedure  Is  capable  of  handling  correlated  Input  random 
variables  (Reference  73)  and  nonsymmetrical  probability  distribution  functions 
(Reference  74),  If  such  parameters  are  known  In  a  problem. 

For  the  purpose  of  the  present  Investigation,  It  will  be  assumed  that  the 
Input  random  variables  are  uncorrelated  and  their  probability  distribution 
functions  are  symmetrical.  The  assumption  Is  motivated  by  the  fact  that  the 
exact  nature  of  these  parameters  Is  seldom  known.  Therefore,  In  the  subse¬ 
quent  analysis.  Equations  3.27-3*30  will  be  used  for  probabilistic  wave  propa¬ 
gation  analyses. 

With  E(Y)  and  V(Y)  known,  the  value  of  the  function  at  one  standard 
deviation  above  (Y+)  and  below  (Y-)  Its  mean  value  can  be  determined  from  the 
following  relation: 

Y*  -  E(Y )  ±  [V(Y)]y*  (3.3D 

We  can  now  proceed  to  apply  the  Rosenblueth  procedure  to  the  wave  propa¬ 
gation  theory  developed  In  the  previous  sections  in  order  to  account  for  the 

randomness  in  the  parameters  £  and  v  .  In  this  connection,  we  will  denote 

£1 

£  and  v  as  the  expected  values  of  these  variables  and  o.  and  a  as 
3  Jt  V 

a 

their  standard  deviations.  The  values  of  these  variables  at  one  standard 
deviation  above  (P)  and  below  (M)  their  mean  values  then  become 


lP  -  l  ♦  o„ 


„M  7 
l  »  £  -  o„ 


P  -  x 
v,  «  vo  +  0 
a  a  v 


M  - 
=  v„  -  c 
a  a  v. 


(3*32) 


Four  deterministic  wave  propagation  calculations  are  conducted  for  the 

,P  .M 


P  M 

v  and  v  ,  as  stipulated  by 

di  3 


four  possible  combinations  of  £  ,  £ 

Equation  3.30.  The  output  from  these  calculations  is  combined  at  successive 
times  (at  a  selected  depth)  according  to  Equations  3*27-3*29  to  calculate  the 
expected  value  and  the  variance  of  each  of  the  dependent  variables  (accelera¬ 
tion,  velocity,  stress,  etc.).  Equation  3*31  can  then  be  used  to  construct 


the  time  histories  of  the  expected  value  and  the  one-standard-devlatlon  bounds 
of  these  variables.  The  computer  code  has  been  developed  to  also  allow  such 
probabilistic  calculations  to  be  made.  To  process  the  results,  the  program 
first  calculates  an  expected  value  for  the  arrival  time  of  the  wave  at  any 
selected  depth  using  the  arrival  time  data  from  the  four  Individual  determin¬ 
istic  calculations.  The  program  then  translates  (shifts)  all  the  waveforms  to 
this  common  arrival  time  for  processing.  Figures  3.22-3.24  Illustrate  some 
typical  probabilistic  results  for  the  acceleration,  stress  and  velocity.  Each 
figure  shows  the  expected  value  and  the  one-standard-devlatlon  bounds  for  each 
wave  form. 
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Figure  3.2.  Typical  two-dimensional  Voronoi  cells;  Reference  77. 
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Figure  3.8.  Amplitude  attenuation  versus  distance  for  periodic  volume 
distribution. 
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Figure  3.10.  Average  wave  speed  versus  distance  for  exponential  volume 
distribution. 
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Figure  3.11.  Amplitude  attenuation  versus  distance  for  exponential  volume 
distribution. 
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Figure  3.12.  Actual  wave  speed  versus  distance  for  combined  periodic- 
exponential  volume  distribution. 
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Figure  3.15.  Wave  profile  construction. 
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Figure  3.16.  Four  wave  profile  example. 
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Figure  3.17.  Postulated  variation  of  material  parameter  as 

of  stress. 
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Figure  3.18.  Particle  acceleration  profile  for  a  periodic  volume 


distribution  (Material  PI)  at  X  -  0.125  in. 
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Figure  3.19.  Stress  profile  for  a  periodic  volume  distribution 
(Material  PI)  at  X  a  0.125  in. 
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Figure  3.21.  Particle  displacement  profile  for  a  periodic  volume 
distribution  (Material  PI)  at  X  »  0.125  in. 
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Figure  3.22.  Probabilistic  results  for  the  particle  acceleration  profile  at 
X  «  0.125  in;  mean  response_with  its  one-standard-deviation 
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Figure  3.23  Probabilistic  results  for  the  stress  profile  at  X  =  0.125  in  ; 

mean  response  with  its  one-standard-deviation  bounds.  Input 
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Figure  3.24.  Probabilistic  results  for  the  velocity  profile  at  X  -  0.125  in  ; 

mean  response  with  its  one-standard-deviation  bounds.  Input 
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CHAPTER  4 


PARAMETRIC  STUDIES 


4.1  GENERAL 

This  chapter  will  present  a  variety  of  typical  results  from  the  developed 
computer  code,  MIC1D.  Depending  upon  which  volume  distribution  function  Is 
used  In  the  modeling,  a  wide  variety  of  predicted  output  can  result  by  varying 
particular  constitutive,  mlcrostructural ,  and  Input  parameters.  The  constitu¬ 
tive  parameters  In  the  wave  propagation  theory  are  ,  aee  ,  A£  ,  and 

A  .  The  mlcrostructural  parameters  are  v  ,  v.  ,  i  ,  B  ,  Y  , 
ee  a  b  o 

Y  ,  and  M  .  The  input  parameters  are  aQ  and  AtQ  .  For  the  purpose  of  the 
parametric  calculations,  the  numerical  values  of  the  constitutive  parameters 
were  kept . constant  for  each  of  the  three  volume  distribution  functions  that 
were  used.  The  parametric  calculations  concentrated  only  on  the  variation  In 
the  mlcrostructural  and  Input  parameters.  In  principle,  however,  the  consti¬ 
tutive  parameters  would  be  a  function  of  the  mlcrostructural  parameters  (see 
Equation  2.6).  Explicit  relationships  for  these  parameters  have  not  yet  been 
determined.  Space  limitations  in  this  report  preclude  presenting  many  cases; 
consequently,  only  major  model  parameters  will  be  considered.  The  chapter  is 
divided  into  three  major  sections  dealing  with  (1)  depth  dependent  behavior, 
(2)  wave  profiles,  and  (3)  probabilistic  profiles.  Subsections  within  each  of 
these  sections  then  present  specific  effects  of  parametric  variation. 

4.2  DEPTH  DEPENDENT  BEHAVIOR 

This  section  will  present  the  effects  of  the  mlcrostructural  and  input 
parameters  on  the  variation  of  the  average  wave  speed  and  amplitude  attenua¬ 
tion  with  depth.  These  wave  propagatlonal  characteristics  are  for  the  case  of 
a  single  wave  moving  Into  regions  which  are  at  rest  in  their  reference  config¬ 
uration. 

4.2.1  Periodic  Volume  Distribution  Case 

For  the  periodic  distribution  model,  the  mlcrostructural  parameters  are 
the  average  porosity  v  and  the  grain  or  characteristic  length  l  .  The 
basic  Input  parameter  Is  the  Initial  amplitude  aQ  .  a  typical  volume  distri¬ 
bution  plot  for  this  case  is  shown  in  Figure  3.1.  Figures  4. 1-4.6  Illustrate 
typical  code  output  for  a  variety  of  parametric  variations.  Figure  4.1 


presents  the  variation  of  average  wave  speed  with  v  .  As  expected,  the 

average  wave  speed  decreases  with  Increasing  porosity  (l.e.,  decreasing  v  ). 

Figure  4.2  shows  that  the  wave  speed  will  Increase  as  l  Increases.  This 

result  Is  apparently  related  to  the  fact  that,  with  an  Increase  In  Z  ,  the 

wave  will  see  fewer  mtcrostructural  changes  per  unit  length  of  travel  and, 

hence,  less  dispersion.  Figure  4.3  Illustrates  the  amplitude  behavior  for 

s 

three  different  Initial  amplitudes,  a  ■  5  x  10  , 

*  4  2  0 

1  x  10  ,  and  5x10  In/s  .  Clearly,  the  expected  result  can  be  seen  in 

that  higher  Initial  amplitudes  decay  faster  than  the  lower  amplitude  waves. 
Figure  4.3,  In  conjunction  with  Figures  4.4  and  4.5,  portray  the  effect  of 
v  on  amplitude  attenuation.  It  Is  observed  that  the  attenuation  rate  Is 

3 

strongly  dependent  on  v  .  As  v  decreases  (l.e.,  Increasing  porosity)  the 
rate  of  attenuation  Increases.  This  result  Is  also  consistent  with  the  varia¬ 
tion  In  wave  speed  with  v  shown  In  Figure  4.1.  Finally,  Figures  4.6  and 
4.4  demonstrate  the  effect  of  Z  on  amplitude  attenuation.  These  figures 
Indicate  that  larger  values  of  Z  result  In  less  attenuation,  which  Is 
consistent  with  the  previous  observation  regarding  the  variation  of  wave  speed 
with  Z  . 

4.2.2  Exponential  Volume  Distribution  Case 

The  exponential  volume  distribution  model  contains  the  mlcrostructural 
parameters  of  the  free  surface  porosity  and  the  depth  rate  of  consolida¬ 

tion  B  .  As  before,  the  Input  parameter  Is  the  Initial  amplitude 
aQ  .  Figures  4.7-4.10  show  typical  results  concerning  the  effects  Qf  these 
parameters  on  the  wave  propagation  variables.  Figure  4.7  shows  the  variation 
of  the  volume  distribution  function  with  distance  for  two  different  values  of 
vb  .  Figure  4.8  shows  the  effect  of  on  the  average  wave  speed.  For  this 

case,  the  wave  speed  Increases  with  depth  due  to  the  overall  decrease  In 
porosity  with  depth.  It  Is  also  apparent  that  an  Increase  In  porosity 
produces  a  slower  wave  speed.  Figures  4.9  and  4.10  show  the  effect  of  vb 
and  the  Initial  amplitude  on  wave  attenuation.  These  results  give  trends 

similar  to  the  previous  observations  for  the  periodic  distribution  function. 

* 

That  Is,  higher  Initial  amplitude  waves  attenuate  faster  and  the  attenuation 
rate  Increases  with  porosity. 

4.2.3  Periodic-Exponential  Volume  Distribution  Case 

For  the  combined  periodic-exponential  distribution  model,  all  four  micro- 


structural  parameters  v  ,  .  i  ,  and  B  are  present.  This,  along  with 

the  Initial  amplitude  aQ  ,  provides  considerable  parameter  variations.  Only 

a  portion  of  the  possible  parametric  variations  will  be  presented,  and  these 

are  shown  In  Figures  4.11-4.14.  Figure  4.11  Illustrates  the  variation  of 

volume  distribution  function  with  distance  for  this  model  for  two  different 

* 

values  of  .  This  combined  function  has  both  oscillatory  and  monotonlc 
depth- dependent  features.  Figure  4.12  shows  the  effect  of  on  the  average 

wave  speed.  Figures  4.13  and  4.14  show  the  effect  of  the  Initial  amplitude 
and  on  wave  attentuatlon.  These  results  portray  the  same  trends  as 

observed  In  the  previous  sections. 

4.3  WAVE  PROFILES 

Time  profiles  of  the  particle  acceleration,  velocity,  and  displacement, 

along  with  the  stress  at  selected  depths  Into  the  medium,  are  presented  In 

this  section.  The  profile  construction  procedure  was  discussed  earlier  In 

Section  3.4.  Only  results  from  the  periodic  volume  distribution  case 

(corresponding  to  model  Material  PI)  will  be  presented:  however,  the  other 

volume  distributions  will  produce  similar  results.  Referring  to  the  wave 

coupling  aspects  discussed  In  Section  3 - ^ *  see  Equation  3.25,  this  section 

will  present  both  uncoupled  (M  -  0)  and  coupled  results  (M  *  0).  The  Input 

_e 

acceleration  profile  used  equal  time  spacing  of  AtQ  =4X10  s  for  all 
cases  presented  here. 

4.3.1  Uncoupled  Results 

Figures  4.15-4.18  Illustrate  uncoupled  results  for  the  four  profiles  at 
two  different  depths,  X  =  0.0  In  ,  X  »  2.5  In  ,  and  X  =  6  In  .  The  Input 
acceleration  wave  Is  shown  In  Figure  4.15  corresponding  to  X  =  0.0  In  .  It 
should  be  pointed  out  that,  for  this  case,  all  Individual  waves  In  a  given 
profile  propagate  Independently  of  each  other.  However,  as  discussed  earlier, 
In  a  given  profile  higher  amplitude  waves  attenuate  faster  than  the  lower 
amplitude  waves. 

4.3.2  Coupled  Results 

For  the  coupled  case.  Equation  3.25  Is  In  effect  and  the  parameter  M 

plays  a  significant  role  In  determining  the  amount  of  coupling.  The  wave 

profile  results  for  the  coupled  case  are  shown  In  Figures  4.19-4.22  for  a 

2 

value  of  M  -  0.04  In  /lb  and  with  v  =  0.85  .  Coupling  effects  through 

-  a  _ 


variation  of  v  (see  Equation  3.25)  for  the  periodic  distribution  case  will 
produce  less  attenuation  than  the  corresponding  uncoupled  results.  This 
occurs  since  Increases  in  v  produce  a  material  with  less  average  porosity 

3 

and,  hence,  dispersive  effects  will  be  reduced. 

*1.4  PROBABILISTIC  PROFILES 

This  final  section  shows  results  of  some  typical  probabilistic  computer 
runs.  The  theoretical  development  was  discussed  previously  in  Section  3.5. 
The  probabilistic  results  consisting  of  the  expected  value  and  the  one- 
standard-deviation  bounds  of  particle  motion  and  stress  are  shown  in 
Figures  4.23-4.26  for  the  case  of  zero  coupling.  These  results  correspond  to 
model  Material  PEI  (the  combined  periodic-exponential  volume  distribution 
case).  Only  the  parameter  A  was  considered  to  be  random  for  these 
calculations.  Therefore,  the  dispersion  in  the  output  quantities  is  only  due 
to  uncertainties  in  l  . 
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Figure  4.1.  The  effect  of  v  on  the  average  wave  speed  versus  depth  for 
a  periodic  volume  distribution  (Material  Fl  with 
l  -  0.1  in). 
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The  effect  of  t  on  the  average  wave  speed  versus  depth 
for  a  periodic  volume  distribution  (Material  PI  with 


Amplitude  Ratio 


1 .0  • 
0.9  - 
O.S  - 

l 

0.7  «> 

i 

0.6  •!, 


0.5 


1.  a  *  5  x  10''  in/s“ 

2.  a  *  1  x  104  in/s 

o  I  ^ 

3.  a  *  5  x  l(r  in/s 

o 


0.4 

0.3 


0.2 


0.0 


0.0 


2.0 


4.0 


6.0 


Propagation  Depth  (in) 


Figure  4.3.  Amplitude  ratio  versus  depth  for  a  periodic  volume 
distribution  (Material  PI  with  l  -  0.1  in  and  v 
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Figure  4.6.  Amplitude  ratio  versus  depth  for  a  periodic  volume 
distribution  (Material  PI  with  t  ■  0.2  in  and  v. 
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Figure  4.8.  Average  wave  speed  versus  depth  for  an  exponential  volume 
distribution  (Material  El  with  B  -  10  in2/lb). 
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Figure  4.9.  Amplitude  ratio  versus  depth  for  an  exponential  volume 
distribution  (Material  El  with  B  ■  10  in^/lb  and 
vfa  -  0.65).  t 
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Figure  A. 10.  Amplitude  ratio  versus  depth  for  an  exponential  volume 
distribution  (Material  El  with  B  ■  10  in2/lb  and 
vfe  -  0.75). 
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Figure  4.12. 


Average  wave  speed  versus  depth  for  a  perodic-exponential 
voloroe  distribution  (Material  PEI  with  t  ■  0.1  in). 


66 


Amplitude  Ratio 


! 


i 


Figure 


* 


Propagation  Depth  (in) 


4.13.  Amplitude  ratio  versus  depth  for  a  periodic-exponential 
volume  distribution  (Material  PEI  with  l  m  0.1  in  and 
v.  -  0.65). 
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Figure  4.14.  Amplitude  ratio  versus  depth  for  a  periodic-exponential 
volume  distribution  (Material  PEI  with  t  *  0.1  in  and 
vfe  -  0.75). 
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Figure  4.15.  Uncoupled  particle  acceleration  profiles  at  various  depths 
(Material  PI  with  v  -  0.85  and  l  -  0.1  in). 
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Figure  A. 17. 


Uncoupled  particle  displacement  profiles  at  various  depths 

(Material  PI  with  v  »0.85  and  •£  *  0.1  in). 
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Figure  4.21.  Coupled  particle  displacement  profiles  at  various  depths 
(Material  PI  with  v  -  0.85  ,  l  -  0.1  in  ,  and 
M  -  0.04  in2/lb) .  o 
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Figure  4.22. 


Coupled  stress  profiles  at  various  depths  (Material  PI  with 
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Figure  4.23.  Probabilistic  acceleration  profiles  at  various  depths;  mean 
response  with  its  one- standard-deviation  bounds  for 
Z  -  0.1  in  and  ■  0.03  in. 
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Figure  4.25.  Probabilistic  displacement  profiles  at  various  depths;  mean 
response  with  its  one-standard-deviation  bounds  for 
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CHAPTER  5 


SUMMARY  AND  RECOMMENDATIONS 


5.1  SUMMARY 

A  one- dimensional  computer  program  .-(referred  to  as  MIC1D)  has  been 
developed  for  analyses  of  explosive  wave  propagation  in  granular  materials 
with  microstructure.  The  theoretical  foundation  of  the  computer  program  is 
based  on  the  distributed  body  concept  advanced  by  Goodman  and  Cowin 
(Reference  31 )  and  the  associated  wave  propagation  studies  conducted  by 
Nunziato,  Walsh,  et  al.  (References  57-63).  The  computer  program  allows  for 
(1)  arbitrary  surface  airblast  loading,  (2)  depth-dependent  volume  distribu¬ 
tion  function  simulating  gravity  effects  in  a  granular  mass,  and  (3)  treatment 
of  grain  size  and  local  porosity  as  random  variables.  Three  forms  of  depth- 
dependent  volume  distribution  functions  are  Incorporated  in  the  program,  l.e., 
a  periodic  form,  an  exponential  form,  and  a  combined  per iodic- exponential 
formulation.  The  user  can  select  any  of  these  forms  for  the  particular 
application  at  hand.  Probabilistic  treatment  of  grain  size  and  local  porosity 
is  accomplished  by  using  a  moment- generating  procedure  due  to  Rosenblueth 
(Reference  73).  The  computer  program  calculates  the  expected  value  and  the 
variance  of  the  output  quantities,  such  as  stress,  particle  motion,  etc.,  due 
to  the  randomness  In  these  variables. 

Application  of  the  computer  program  is  demonstrated  by  presenting  the 
results  of  a  series  of  parametric  calculations  dealing  with  propagation  of 
acceleration  waves  In  granular  media.  It  is  shown  that,  within  the  range  of 
variables  studied,  local  porosity  and  grain  size  characteristics  (which  are 
reflected  In  the  mlcrostruct.ural  parameters  va  ,  vfe  ,  B  ,  l  ,  fQ  ,  Y) 
play  an  Important  role  In  wave  attenuation  in  such  materials.  The  effect  of 
grain  size  parameter  i  on  wave  attenuation  during  short  propagation 
distances  (less  than  a  hundred  grains)  is  Interesting  and  requires  experimen¬ 
tal  verification. 

5.2  RECOMMENDATIONS 

Several  recommendations  are  made  for  further  development  In  the  analyti¬ 
cal  aspects  and  for  the  experimental  validation  of  the  present  work.  First, 
the  wave  propagation  model  should  be  validated  with  an  experimental  setup  that 
captures  the  dynamic  load  transfer  and  wave  motion  across  several  grain 


boundaries.  This  can  be  accomplished  via  dynamic  photoelasticity  as  described 
by  Rossmanlth  and  Shukla  (Reference  75).  For  example,  the  dynamic  photoelas- 
tlc  laboratory  at  the  University  of  Rhode  Island  can  take  data  at  a  rate  of 
10^  frames/sec.  It  Is  possible  to  measure  wave  speed  and  amplitude  attenua¬ 
tion  as  the  dynamic  loads  are  transferred  across  a  finite  number  of  particles. 
Second,  the  coupling  effects  for  constructing  wave  profiles  from  a  group  of 
traveling  waves  should  be  examined  further.  In  the  present  formulation,  the 
coupling  effects  due  to  spatial  redistribution  of  local  porosity  as  the  waves 
propagate  through  the  medium  were  neglected.  The  consequence  of  this 
assumption  needs  to  be  examined.  Third,  in  the  present  study  the  constitutive 
moduli  were  taken  to  be  constant,  but,  In  actuality,  they  depend  on  the  basic 
ml cr ©structural  parameters.  Attempts  should  be  made  to  relate  these  moduli  to 
some  basic  mlcrostructural  variables  of  the  granular  medium.  Studies  In 
fabric  tensor  models  proposed  by  Oda,  Nemat-Nasser ,  et  al . ,  (Reference  67)  may 
provide  some  direction  for  developing  the  desired  relationships.  Fourth,  the 
probabilistic  analysis  In  this  study  was  only  relevant  to  the  periodic  aspects 
of  the  volume  distribution  model  since  only 

v  and  t  were  treated  as  random  variables.  The  present  probabilistic 
d 

analysis  should  be  extended  to  treat  in  the  exponential  distribution 

model  also  as  a  random  variable.  Finally,  the  wave  propagation  model  should 
be  extended  to  two-dimensional  geometry.  The  variation  In  local  porosity  In 
two  dimensions  can  then  be  Incorporated  In  such  a  formulation  without  recourse 
to  a  probabilistic  analysis. 
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APPENDIX  A 

WAVE  PROPAGATION  COMPUTER  CODE 


APPENDIX  A 


GLOSSARY 
NAME 
A(N, J) 

A 1  (  N  ,  J  ) 

A2 ( N , J ) 
A3(N, J) 
A4(N, J) 

AB ( 2 J ) 

AFF 
AFFF 
AM  (  2  J  ) 
AP(2J; 
APLOT ( 2 J 
AR(N, J) 

BB 

C1(N,J) 
C2(  N  ,  J  ) 
C3(N, J) 
DD(J,  .) 

DELTR 

DELX ( N , J ) 
DISP ( N , J ) 

DM(  J ) 


WAVE  PROPAGATION 
COMPUTER  CODE 


OF  MAJOR  VARIABLE  NAMES 


DEFINITION 


Accelerat.1  on 

array 

for 

N  depths  and  J 

waves 

Acceleration 

array 

for 

probabilistic 

cal cul at  1 ons 

Acceleration 

array 

for 

pr obab 11 1st! c 

calculations 

Acceleration 

array 

for 

probabil isti c 

calculations 

Acceleration 

array 

for 

probabil istic 

calculations 

Amplitude  array  for  acceleration  plotting 
<xe  >  material  constant 
otee  .  material  parameter 

Acceleration  profile  minus  one-standard-deviation 
Acceleration  profile  plus  one-standard-deviation 
.)  Acceleration  array  for  probabilistic  plotting 
Normalized  amplitude  array 

B,  power  of  exponential  volume  distribution'  function 
y  coefficient  in  amplitude  equation 
k  coefficient  in  amplitude  equation 

Array  for  volume  distribution  function  | 

I 

Particle  displacement  array  for  probabilistic 

plotting  i 

i 

Initial  time  Increment  between  input  acceleration 

waves  I 

I 

Distance  between  two  consecutive  waves.  j 

i 

Array  for  particle  displacement  j 

Particle  displacement  profile  minus  one- standard-  * 

deviation  J 

i 

t 


DP(J)  Particle  displacement  profile  plus  one-standard- 

deviation 

DPLOT(J)  Displacement  array  for  plotting 

D1(N,J)  Displacement  array  for  probabilistic  calculations 

D2(N,J)  Displacement  array  for  probabilistic  calculations 

D3(N,J)  Displacement  array  for  probabilistic  calculations 

D4(N,J)  Displacement  array  for  probabilistic  calculations 

EA(2J)  Expected  acceleration  profile  array 

ED ( J )  Expected  particle  displacement  profile  array 

ENUAA  (va)0  value  of  va  in  the  reference  state 

ENUB  vfc  material  function  for  exponential  volume  distri¬ 

bution  function 


ES  ( J ) 

Expected  stress  profile  array 

EV  ( J ) 

Expected  particle  velocity  profile  array 

EX 

M,  material  constant  in  equation  (3.25) 

FLFF 

A  ,  material  constant 
e 

FLFFF 

A  ,  material  constant 
ee 

GG 

Y ,  specific  weight  of  the  material  used  in  equation 
(3.2) 

GL 

l,  characteristic  length  for  periodic  distribution 
form  (3.1) 

GO 

Y0 ,  weight  density  of  the  material 

H 

S./40,  Interval  measure  for  numerical  integration 

J 

Looping  index  for  number  of  input  waves 

L 

Maximum  limit  of  the  looping  parameter  N 
gation  depth 

for  propa- 

M 

Looping  index  for  probabilistic  calculations 

N 

Looping  index  for  depth  incrementing 

NR 

Maximum  number  of  runs  for  probabilistic 
calculations.  (NR-1  yields  deterministic 

case ) 

A- 3 

NV 


Value  of  looping  parameter  at  a  desired  depth 


NW  Number  of  Initial  Input  wave  amplitudes  used  to 

construct  a  profile 

PSTRES(J)  Array  used  for  plotting  of  stress  profile 

SIGNUA  ova*  the  standard  deviation  for  volume  distribution 

parameter  va 

SGL  o  * i  the  standard  deviation  for  characteristic  length, 

r 

SM ( J )  Stress  profile  minus  one-standard-deviation 

SP ( J )  Stress  profile  plus  one-standard-deviation 

SPLOT ( J ,  . )  Stress  array  for  probabilistic  plotting 
STRES ( N , J )  Wave  front  stress  array 

S1(N,J)  Stress  array  for  probabilistic  calculations 

S2(N,J)  Stress  array  for  probabilistic  calculations 

S3 ( N , J )  Stress  array  for  probabilistic  calculations 

S4(N,J)  Stress  array  for  probabilistic  calculations 

TAVG(J)  Average  time  for  probabilistic  run 

T(N,J)  Time  of  propagation 

TSUM ( N , J )  Arrival  time  increment 

TP(J)  Time  array  used  for  profile  plotting 

TPLOT(J)  Time  array  used  for  profile  plotting 

T1(N,J)  Time  array  for  probabilistic  calculations 

T2(N,J)  Time  array  for  probabilistic  calculations 

T3(N,J)  Time  array  for  probabilistic  calculations 

T4 ( N , J )  Time  array  for  probabilistic  calculations 

UBAR  Average  wave  speed  as  calculated  in  subroutine 

Gauss 

UF(N)  Fast  wave  speed 

UGRN  Initial  wave  speed  on  the  free  surface  of  the  medium 

at  X-0.  Also  equal  to  granule  wave  speed. 


U(N,J)  Average  wave  speed 

UREAL(N,J)  Actual  wave  speed 

V(N,J)  Particle  velocity 

VM ( J )  Particle  velocity  profile  minus  one-standard- 

deviation 

VP(J)  Particle  velocity  profile  plus  one-standard- 

deviation 

VPART(J)  Particle  velocity  for  profile  plotting 

VV( J , . )  Particle  velocity  array  for  probabilistic  plotting 

VI ( N , J )  Particle  velocity  array  for  probabilistic  calcula¬ 

tions 

V2(N,J)  Particle  velocity  array  for  probabilistic  calcula¬ 

tions 

V3(N,J)  Particle  velocity  array  for  probabilistic  calcula¬ 

tions 

V4 ( N  f  J )  Particle  velocity  array  for  probabilistic  calcula¬ 

tions 

X ( N )  Depth  into  the  media 

YAM ( J )  Initial  input  amplitudes 


‘5  V 


ft; 

S* 

ii 


i 

fit 

V*» 

ft*' 


V  V. 


A- 5 


WAVE  PROPAGATION  CODE 
FLOW  CHART 


Choose  v-Function  From  Menu 

1.  Periodic 

2.  Exponential 

3.  Periodic-Exponential 


Call  Subroutine  RUNKUT  to  Calculate 


EXAMPLE  RUN 


As  an  example  of  the  use  of  the  program  MIC1D,  a  typical  run 
will  now  be  presented.  The  first  prompt  from  the  program  which 
will  appear  on  the  CRT  screen  will  be  a  menu  for  the  type  of 
volume  distribution  to  be  used,  l.e. 

|  CHOOSE  THE  TYPE  OF  VOLUME  DISTRIBUTION  | 

|  FUNCTION  TO  BE  USED  IN  THE  ANALYSIS  FOR  | 

|  THE  GRANULAR  MEDIUM  IN  CONSIDERATION  | 

ENTER  THE  PROPER  SELECTION  NUMBER  PLEASE 


1.  PERIODIC  VOLUME  D ISTR IBUTON  FUNCTION  f 

I 

2.  EXPONENTIAL  VOLUME  DISTRIBUTION  FUNCTION  | 

I 

3.  PERIODIC-EXPONENTIAL  COMBINED  FUNCTION  | 


The  user  should  respond  with  the  appropriate  choice. 

The  name  of  the  Input  data  file  is  requested  next  through 
the  prompt 


ENTER  NAME  OF  THE  DATA  FILE.  ( NAME. DAT ) 
PE OAT. 


The  program  requires  the  Input  data  AFF,  AFFF,  BB,  DELTR,  ENUAA, 
ENUB,  EX,  FLFF ,  FLFFF ,  GG ,  GL ,  GO,  L,  NR,  SIGNUA,  SGL  (see 
glossary  for  description  of  each  of  these  quantities).  The  data 
Is  read  in  an  unformated  fashion,  and  a  typical  data  file  would 
look  like 


~3.E5.5.E8.30.0,4.E-6,0.996,0.7,0.0,3.E5,-750.0,7.2E-2,0. 1 , 
2. 4E-4, 1000, 1,0. 0,0. 03 


A-9 


Depending  upon  which  volume  distribution  that  i3  selected  and 
whether  the  run  is  deterministic  or  not,  some  of  the  input  data 
will  not  actually  be  used. 

Next  the  number  of  Initial  wave  amplitudes  are  requested  by 

INPUT  NUMBER  OF  INITIAL  WAVE  AMPLITUDES  TO  BE  USED 

Finally  the  code  requests  the  values  of  each  of  the  initial 
input  waves 


ENTER  INITIAL  VALUES  OF  THE  WAVE  AMPLITUDES 


The  input  is  now  complete  and  the  code  will  print  out  all  of 
this  data  so  that  the  user  can  check  to  see  if  the  input  has  been 
done  correctly.  A  typical  output  of  this  step  is  shown  below  for 
the  case  of  the  previous  data  file  and  with  three  initial 
ampl i tudes 

INPUT  DATA 

AFF  =  -300000.0 

AFFF  =  5 . 0000000E+08 

BB  =  30.00000 

DELTR  =  4.0000000E-06 

ENUAA  =  0.9960000 

ENUB  =  0.7000000 

EX  =  O.0000000E+00 

FLFF  =  300000.0 

FLFFF  =  -750.0000 

GG  =  7 . 1 999997E-02 

GI.  =  0.1000000 

GO  *  2 . 3999999E-04 

L  =  1000 

AMPLITUDE!  1)  =  0.10000E+04 

AMPLITUDE!  2)  =  0.50000E+04 

AMPLITUDE!  3)  =  0.10000E+05 


PRESS  RETURN  TO  CONTINUE 


Once  the  "return  key"  has  been  depressed,  the  program  will 
start  its  computation  and  the  following  will  typically  be  seen  on 
the  screen 

PROGRAM  16  NOU  RUNNING 

COMPUTATION  OF  DELTA  X  BETWEEN  NAVES  ENDS 

J«  2 

N=  644 

T<N,J>=  7 . 3260060E-05 

MAXIMUM  LIMIT  =  7 . 3219308E-05 

COMPUTATION  OF  DELTA  X  BETWEEN  WAVES  ENDS 

J=  3 

N=  888 

T<N,J>=  7 . 3226467E-05 

MAXIMUM  LIMIT  *  7 . 3219308E-05 

The  message  concerning  the  computation  of  DELTA  X  is  related  to 
the  theory  in  section  3.1*  dealing  with  the  stress  calculations. 
This  computation  will  stop  before  the  completion  of  the  entire 
N-loop,  and  this  comment  lets  the  user  know  when  this  occurs. 

The  above  example  is  for  the  case  of  three  Initial  wave 
amplitudes  with  L-1000. 

When  the  code  is  finished  with  its  calculations,  the  user  is 
prompted  with  the  following  menu  for  output  results 

CHOOSE  THE  TYPE  OF  GRAPHICS  OPTION 

1.  DETERMINISTIC  PROFILE  PLOT 

2.  PROBABILISTIC  PROFILE  PLOT 


3.  DEPTH  DEPENDENT  PLOT 


Selecting  for  example  item  #3,  depth  dependent  plots,  the 
final  plotting  menu  will  appear. 

ENTER  ITEM  NUMBER  FOR  DEPTH  DEPENDENT  PLOT 

1.  ACTUAL  WAVE  SPEED 

2.  AVERAGE  WAVE  SPEED 

3.  AMPLITUDE  BEHAVIOR  WITH  DEPTH 

4.  VOLUME  DISTRIBUTION  FUNCTION 

5.  MU-COEFFICIENT 

6.  KAPPA-COEFFICIENT 

7.  NO  GRAPHICS  i.e.  DATA  OUTPUT  IN  A  FILE 

Specifics  on  the  final  stages  of  plotting  will  not  be  given 
since  these  will  be  system  and  software  dependent  and  thus  will 
vary  from  system  to  system.  An  example  output  corresponding  to 
the  example  input  data  shown  previously,  is  given  in  Figure  A.l. 


Vm  8p««4  (lda/i) 


AVERAGE  WAVE  SPEED 


PrepM*tlen  Depth  (In) 


Figure  A-l.  Plotted  output  for  example  run  (continued). 
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Amplitude 


AMPLITUDE  BEHAVIOR  WITH  BERTH 


Propagation  Depth  (in) 


Figure  A-l.  Plotted  output  for  example  run  (concluded). 
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PROCRAM  M  I  C  1  D 


ONE  DIMENSIONAL  WAVE  PROPAGATION  CODE  FOR  MICROSTRUCTURAL 
MATERIALS  MODELLED  BY  COODMAN-COWIN  DISTRIBUTED  BODY 
THEORY.  CODE  EMPLOYS  SINGULAR  SURFACE  WAVE  THEORY 
ORIGINALLY  DEVELOPED  BY  NUNZIATO  tt  ml. 

WRITTEN  BY 

PROFESSOR  MARTIN  H.  SADD 


MOHAMMAD  N.  HOSSAIN 


DEPARTMENT  OF  MECHANICAL 
ENGINEERING  Si  APPLIED  MECHANICS 
UNIVERSITY  OF  RHODE  ISLAND 
KINOSTON,  RI  02881 
SEPTEMBER.  1986 

WRITTEN  FOR 


DR.  BEHZAD  ROHAN I 

US.  ARMY  ENGINEER 

WATERWAYS  EXPERIMENT  STATION 
VICKSBURG.  MS  39180 


DIMENSION  SIZES  ARE  RELATED  TO  DEPTH  COUNTER  "N"  (1<N<L>  AND 
NUMBER  OF  WAVES  COUNTER  "J”  (1<J<NW>»  MAXIMUM  L-3000  AND  , 
MAXIMUM  J-15 

DIMENSION  A ( 3000.  15). A1 (3000. 15).  A2(3000.  15),  A3(3000.  15),  A4 
»  (3000.  15), AB ( 30 ) . AM ( 30 ) ,  AP ( 30 ) ,  APLOT ( 30.  3),  AR(3000.  15). Cl (3000, 
»  15).  C2 ( 3000.  IS.  C3 ( 3000.  15),  01(3000.  15).  D2(3000,  15),  D3(3000.  15). 
»  D4( 3000,  15). DELX ( 3000.  15),  DD( 15.  3),  DISP (3000.  15), DPLOT< 15). 

*  DP ( 1 5 ) ,  DM (15).  EA ( 30 ) ,  ED (15),  ES (15),  EV (15).  NMAX (15), 

»  PSTRES (15). STRES ( 3000.  15>,S1(3000.  15),  82(3000.  15).  83(3000.  IS). 

»  34(3000,  15),  SM (15),SP(1S)<  SPLOT ( 1 5,  3 ) >  T ( 3000,  15), 

*  TAVG(  15).  TK3000,  15).  T2(3000,  15).  T3(3000,  15),  T40000,  15),  TSUM 

»  (3000.  15).  TP(15),  TPLOT (30),  VIOOOO,  15).  V2  0000,  15),  V3(3000,  15). 

*  V4 ( 3000,  15),  V(3000, 15),  VP(20>.  VPART( 15),  VM( 15),  VV( 15.  3).  UREAL 

*  (3000, 15), U(3000> 15). X ( 3000 ) , YAM (15) 

CHARACTER  ANS*B» DATNAME*8, DATA*8, TITLE*50. XLABEL*50, YLABEL*50. 

»  LINLABEL( 10) #25,  LINE *25,  PANS*5 
INTEGER  JJ. KK.  IANS 

REAL  ALO,  AHI,  AMP,  H,  U,  UL,  UH,  X,  Y,  YAM,  YO 

COMMON  /GRAIN/  AFFO,  AFFF,  BB,  DUM.  ENUB,  FLFFO,  FLFFF,  00,  GL.  00.  TIME 
COMMON  /CONTRL/  P.P1.Q, 01 
P  -  1.  0 
PI  «  0.  0 


V  *,-0  9  4  #,*  « 


I  M  »,* 


01  -0.0 
WRITE <6. 2) 

2  FORMAT! '1  V//) 

PRTYT*. ' - 

PRINT*. ' !  CHOOSE  THE  TYPE  OF  VOLUME  DISTRIBUTION  ! 
PRINT*.  ' !  FUNCTION  TO  BE  USED  IN  THE  ANALYSIS  FOR  I 
PRINT*. ' i  THE  ORANULAR  MEDIUM  IN  CONSIDERATION  t 


PRINT*.  ' - ' 

PRINT*.  '  . 

PRINT*. '  ENTER  THE  PROPER  SELECTION  NUMBER  PLEASE' 
PRINT*.  '  ' 

PRINT*.  ' - ' 

PRINT*,'!  1.  PERIODIC  VOLUME  DISTRIBUTON  FUNCTION  !' 

PRINT*.'!  !' 

PRINT*,  ' !  2.  EXPONENTIAL  VOLUME  DISTRIBUTION  FUNCTION  !' 

PRINT*.  ' !  i  ' 


PRINT*. 'I  3.  PERIODIC-EXPONENTIAL  COMBINED  FUNCTION  !' 

PRINT*.  ' - - ' 

PRINT*.  '  ' 

READ! 9. 9)  KK 
9  FORMAT! 12) 

IF!  !KH.  NE.  1).  AND.  !KK.  NE.  2).  AND.  !KK.  NE.  3)  >  GOTO  1 
IF ! KK.  EQ.  1)  THEN 
Q  -  O.  O 
Q1  =  1.0 
END  IF 

IF!KK.  EQ.  2)  THEN 
P  -  0.  O 
PI  -  1.  O 

•  END  IF 

READ  INPUT  DATA  FROM  THE  EXISTING  DATA  FILE 

WRITE!*. *) 'ENTER  NAME  OF  THE  DATA  FILE.  ! NAME.  DAT) ' 

READ! 9.  10)  DATNAME 
10  FORMAT! 1A8) 

OPEN! UNIT— 2,  FILE-DATNAME.  STATUS- 'OLD ' > 

READ! 2.  *)  AFF.  AFFF.  BB,  DELTR, ENUAA, ENUB. EX. FLFF. FLFFF . GG. GL, GO,  L, 

#  NR.  SENUA. SCL 
15  CLOSE! UN IT-2) 

AMPLITUDE  INPUT  SECTION.  ENTER  THE  !NW)  INITIAL  WAVE  AMPLITUDES. 
FINAL  WAVE  !J-NW)  IS  USED  ONLY  AS  AN  END  MARKER  AND  WILL  NOT 
BE  PLOTTED  IN  PROFILE  RESULTS. 

PRINT*.  '  ' 

PRINT*.  'INPUT  NUMBER  OF  INITIAL  WAVE  AMPLITUDES  TO  BE  USED' 

READ!  9,*)  NW 
PRINT*.  '  ' 

PRINT*.  'ENTER  INITIAL  VALUES  OF  THE  WAVE  AMPLITUDES' 

READ! 5. *)  !YAM! J),  J-l.  NW) 

PRINT*. '  ' 

PRINT*. '  INPUT  DATA' 

PRINT*.  '  ' 

WRITE!6.  *>  'AFF  -  '.AFF 
WRITE! 6. *)  'AFFF  -  ' . AFFF 
WRITE!  6.  *)  'BB  -  '.BB 
WRITE! 6.*)  'DELTR  -', DELTR 
WRITE!*,*)  'ENUAA  -  '.  ENUAA 


WRITE (6/  *)  'ENUB  -  '.ENUB 
WR I TE  ( 6.  *)  'EX  ■  '.  EX 
WRITE  (6.  *)  'FLFF  -  ' . FLFF 
WRITE (6. •)  'FLFFF  -  ' »  FLFFF 
WRITE (6.  *>  'GO  -  '.  OG 
WRITE (6.  •)  'GL  -  '.  OL 
WRITE<6»  *)  'GO  >  '.GO 
WRITE (6.  *>  'L  ■  '.  L 
PRINT*.  '  ' 

DO  I-l.NW 

WRITEI6. 800)  I.YAM(I) 

END  DO 
PRINT*.  '  ' 

PRINT*.  'PRESS  RETURN  TO  CONTINUE.  .  . 
READ  <5,  *) 

PRINT*.  '  ' 

PRINT*.  'PROGRAM  IS  NOW  RUNNING' 

AFFO-AFF/CO 

FLFF0*FLFF/C0 

UGRN  *  SORT ( FLFFO ) 

TAVG<  1  )=0.  0 
H-CL/40.  O 

PROBABILISTIC  LOOP  STARTS 

IF  (NR.  EQ.  1 )  GO  TO  17 
PRINT*.  '  ' 

PRINT*. 'PROBABILISTIC  RUN  WITH  DATA' 
PRINT*.  '  ' 

WRITE <6. *)  ' ENUAA  -  ' . ENUAA 
WRITE (6. *)  'EX  = ' , EX  . 

WRITE  (  6.  •)  'NR  «'.  NR 
WRITE (6.  *>  'SIONUA  ■  '.SENUA 
WRITE(6. *)  'SOL  “ ' ,  SOL 
PRINT*.  '  ' 

PRINT*.  'PRESS  RETURN  TO  CONTINUE.  .  .  ' 
READ  (9.  *) 

ENUAP  »ENUAA+SENUA 
ENUAM-ENUAA-SENUA 
CLP  *  GL  *  SGL 
CLM  -  OL  -  SOL 
DO  71  Hal. NR 
IF  (NR.  EQ.  1)  GO  TO  18 
IF  (M  EQ.  1)  THEN 
ENUAA-ENUAP 
GL  «  CLP 
END  IF 

IF  (M.  EQ.  2)  THEN 
ENUAA-ENUAP 
OL  «  CLM 
END  IF 

IF  (M.  EQ.  3)  THEN 
ENUAA  •  ENUAM 
OL  >  GLP 
END  IF 

IF  (M.  EQ.  4)  THEN 
ENUAA  «  ENUAM 
OL  «  OLM 
END  IF 

18  ENUA  -  ENUAA 


DO  70  J  -  l.NU 
DUM-0  0 
Y0»YAM(J> 

DO  69  N  -  If  L 
X(N)«N*CL/40. 

Y-X(N> 

PERFORM  NUMERICAL  INTEGRATION  FOR  THE  TRAN8IT  TIME 
AND  CALCULATE  THE  AVERAGE  WAVE  SPEED 

XL»INT(<N-l>/20.  )*GL/2.  O 
NN-INT< (L-20.  >/20.  ) 

DO  20  1*1. NN 
LL-20.  *1  +  1 

IF  (N.  EG.  LL)  GOTO  30 
20  CONTINUE 
GOTO  40 
30  DUM-TIME 

40  IF  (N.EQ.  1)  TSUM(l.J)  *  H/UORN 

IF  <N.  GT.  1 )  TSUM<N,J>  -  TSUM(N-1.  J)+H/<U(N-1.  J>*1000.  > 

T<N, J)  -  (J-l >*DELTR+TBUM(N»  J) 

IF  (J.  EG.  1)  GO  TO  90 
IF  (T(N,  J>.  CT.  T(L-1.  1)  >  THEN 
PRINT*. '  ' 

PRINT*. 'CAN  NO  LONGER  COMPUTE  DELTA  X  BETWEEN  WAVES' 

PRINT*.  'J-'.  J 
PRINT*.  'N-'.  N 
PRINT*.  'T(N.  J)-'.T(N.  J) 

PRINT*.  'MAXIMUM  LIMIT  -  '.T(L-1,1> 

CO  TO  70 

END  IF  « 

DO  NC*N.  N+200 
NT-NC 

IF  (T(N. J).  LE  T(NC. J-l) )  00  TO  49 
END  DO 

PRINT*.  ' M= ' .  M 
PRINT*.  'N-'.  N 
PRINT*.  'J-'.  J 

PRINT*, 'ERROR  IN  THE  AMPLITUDE  TIME-TEST' 

GO  TO  630 

49  DELX<N,  J-l )«<NT-1 )*H+<T<N,  J)-T<NT-1. J-l ) )*<H/(T(NT.  J-l ) 

#— T<NT-1, J-l)) ) — N*H 
NMAX < J)  -  NT 
MAXN-NMAX( J) 

ASUM-0  0 
DO  JK-l.J-1 

ASUM«ASUM*A ( N.  JK ) 

END  DO 

SI0MA«C3<N.  1 ) *00* ASUM*DEL X ( N.  J-l) 

STRES(N,  J)-SIGMA*STRES<N. J-l ) 

ENUA  -  1,  -  ( 1.  -ENUAA) *EXP<-EX*STRES(N. J) ) 

90  CALL  GAUSS  (ENUA.  J.  KK.  N.  XL.  Y.  UBAR ) 

U(N.  J)-UBAR/1000 

SOLUTION  OF  THE  WAVE  AMPLITUDE  BY  FOURTH  ORDER  RUNOE-KUTTA  METHOD 

CALL  RUNKUTtAMP. ENUA,  H.  J,  N,  PKAPA,  PMU,  VDFUN,  RWSPD,  Y.  YO) 

A(N. J)  -  AMP 

AR(N. J)  -  A(N.  J)/YAM< J) 
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Cl  <N. J)  -  PMU 
C2(N, J)  -  PKAPA 
C3(N. J)  -  VDFUN 
UREAL (N, J)  -  RWSPD 
ASUM-ASUM+A(N,  J) 

IF  <J.  EQ.  1)  THEN 
STRES(N,  J)=0.  0 
V(N.  J>-0.  O 
DISP  ( N»  J>«0.  0 
CO  TO  55 
END  IF 

V(N.  J)=V(N» J-l >+A(N,  J-l )*<T(N.  J)-T<N. J-l ) ) 

DISP (N. J)*  5*( V(N, J)+V(N,  J-l >  >*(T(N,  J)-T(N»  J-l ) >+DISP(N.  J-l) 
IF  (NR.  EQ.  1 )  00  TO  65 
IF  <M.  EQ.  1)  THEN 
A1 (N, J)  ■  A(N>  J) 

T1 (N. J)  -  T(N« J> 

SI <N. J)*STRES(N, J) 

VI (N. J)*V(N,  J> 

D1 (N,  J)«DISP(N,  J) 

END  IF 

IF  (II.  EQ.  2)  THEN 
A2(N, J)  ■  A(N.  J) 

T2(N, J)  -  T(N, J> 

S2(N.  J)-STRES(N.  J) 

V2(N,  J ) »V ( N»  J) 

D2(N, J)»DISP(N.  J) 

END  IF 

IF  (M.  EQ  3)  THEN 
A3(N» J)  -  A(N» J) 

*  T3SN,  J)  ■  T(N,  J) 

S3(N.  J)-STRES(N,  J) 

V3(N. J)-V(N,  J> 

D3(N# J)*DISP(N.  J) 

END  IF 

IF  (M.  EQ.  4)  THEN 
A4(N. J)  -  A(N.  J) 

T4(N. J)  -  T(N.  J) 

S4(N. J)-STRES(N.  J> 

V4(N, J)-V(N.  J) 

D4(N. J>-DISP(N.  J) 

END  IF 
CONTINUE 
CONTINUE 
CONTINUE 
PRINT*. '  ' 

START  OF  FRIENDLY  INTERACTIVE  OPTIONS 

PRINT*, 'CHOOSE  THE  TYPE  OF  ORAPHICS  OPTION' 

PRINT*. '  ' 

PRINT*. '1  DETERMINISTIC  PROFILE  PLOT' 

PRINT*, '  ' 

PRINT*, '2.  PROBABILISTIC  PROFILE  PLOT' 

PRINT*. '  ' 

PRINT*.  '3.  DEPTH  DEPENDENT  PLOT' 

PRINT*, '  ' 

READ (5, *)  IANS 

IF  (  IANS  EQ.  1  )  CO  TO  75 

IF  (  IANS.  EQ.  2  )  CO  TO  135 


DEPTH  DEPENDENT  PLOT ‘ 


CO  TO  75 
CO  TO  135 
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IF  (  IANS.  EO.  3  )  00  TO  80 

IF  ((  IANS.  NE.  1  ).  AND.  (IANS.  NE.  2).  AND.  (IANS.  NE.  3)  >  00  TO  72 
75  PRINT  760 

PRINT  600.  MAXN 
READ (5. 650)  NV 
XCURNT =NV*GL/40 
PRINT  680.  XCURNT 

TIME  PROFILE  PLOTTING  OPTIONS 

PRINT*.  '  ' 

PRINT*. 'ENTER  OPTION  NUMBER  FOR  DETERMINISTIC  PROFILE  PLOT' 
PRINT*.  '  ' 

PRINT*. '1.  PARTICLE-DISPLACEMENT  PROFILE' 

PRINT*.  '  ' 

PRINT*. '2.  PARTICLE-VELOCITY  PROFILE' 

PRINT*.  '  ' 

PRINT*. '3.  PARTICLE-ACCELERATION  PROFILE' 

PRINT*. '  ' 

PRINT*. '4.  STRESS  PROFILE' 

PRINT*.  '  ' 

READ (5. 580)  JF 

CO  TO  (112.115.120,125)  JF 

DEPTH-DEPENDENT  PLOTTING  OPTIONS 

80  PRINT*, 'ENTER  ITEM  NUMBER  FOR  DEPTH  DEPENDENT  PLOT' 

PRINT*. '  ' 

PRINT*. '1.  ACTUAL  WAVE  SPEED' 

PRINT*.  '  ' 

PRINT*. '2.  AVERAGE  WAVE  SPEED' 

PRINT*.  '  ' 

PRINT*. '3.  AMPLITUDE  BEHAVIOR  WITH  DEPTH' 

PRINT**  '  * 

PRINT*. '4.  VOLUME  DISTRIBUTION  FUNCTION' 

PRINT*. '  ' 

PRINT*. '5.  MU-COEFFICIENT' 

PRINT*.  '  ' 

PRINT*. '6.  KAPPA-COEFFICIENT' 

PRINT*. '  ' 

PRINT*.  '7.  NO  GRAPHICS  i.  e.  DATA  OUTPUT  IN  A  FILE' 

READ (5. 580)  JJ 

GO  TO  <100.  105.  110,  130.  140,  150.200)  JJ 

CALL  GRAPHICS  SUBROUTINES  TO  DRAW  THE  COMPUTED  OUTPUT 

PLOT  THE  ACTUAL  WAVE  SPEED 

100  UMIN=200. 

UMAX*0  0 
DO  J-l.NW 
DO  1=1. L 

IF  (UREAL< I. J).  LE.  UMIN)  UMIN=UREAL( I. J) 

IF  ( UREAL ( I. J).  GE.  UMAX)  UMAX=UREAL ( I . J ) 

END  DO 
END  DO 

WRITE(6. 900)  UMIN, UMAX 

PRINT*. 'ENTER  LOWER  &  UPPER  ORDINATES  OF  WAVE  SPEED' 

READ  <  5.*)  VL.VH 

TITLE* 'GACTUAL  WAVE  SPEED*' 


r>  n  n  n  r>  n 


XLABEL* '•Propagation  Depth(in)*' 

YLABEL- '*Wave  Speed  <kin/t>*' 

CALL  WAVGRAF( X,  0  0.  X  (L  ) »  UREAL.  VL.  VH,  L.  1.  3000.  19.  LINLABEL. 
RXLABEL, YLABEL.  TITLE.  .  TRUE.  .  .  TRUE.  ) 

GO  TO  500 

PLOT  THE  AVERAGE  WAVE  SPEED 

109  UMAX-0  0 
DO  J=l, NW 

L I NLABEL (J)= '  ' 

DO  1  =  1.  L 

IF  (U(I.  J)  .  GE.  UMAX)  UMAX=U< I.  J) 

END  DO 
END  DO 

WRITE (6,  1000)  UMAX 
VL=0.  0 

PRINT*.  'ENTER  UPPER  ORDINATE  OF  AVERAGE  WAVE  SPEED' 

READ  (5,*)  VH 
J=1 

DO  WHILE  <  J  .  NE.  0  ) 

PRINT*. 'ENTER  THE  LINE  «  TO  BE  LABELED' 

PRINT*. '(TYPE  0  IF  NO  LABELS  NEEDED  >' 

READ<  5,  *)  J 
IF  (  J  .  NE.  O  )  THEN 
PRINT  640,  J 
READ (5,  570)  LINE 
L I NLABEL ( J ) =L I NE 
END  IF 
END  DO 

TITLE="*AVERACE  WAVE  SPEED*' 

XLABEL= '^Propagation  Depth  (in)*' 

YLABEL* '*Wave  Speed  (kin/*)*' 

CALL  WAVGRAF ( X .  0  0,  X(L).  U.  VL.  VH,  L.  1. 3000.  15.  LINLABEL. 
4XLABEL. YLABEL. TITLE.  .  TRUE.  .  .  TRUE.  ) 

CO  TO  500 

PLOT  THE  AMPLITUDE  BEHAVIOR  WITH  DEPTH 

110  T I TLE=  ' *AMPL I TUDE  BEHAVIOR  WITH  DEPTH*' 

XLABEL* '*Propagat ion  Depth  (in)*' 

YLABEL*  '*Amplitu<ie  Ratio*' 

DO  J*l, NW 

LINLABEL (J)  = '  ' 

END  DO 
J-l 

DO  WHILE  (  J  .  NE.  O  ) 

PRINT*,  'ENTER  THE  LINE  #  TO  BE  LABELED' 

PRINT*. '(TYPE  0  IF  NO  LABELS  NEEDED  >' 

READ (5,  *)  J 
IF  (  J  .  NE.  O  )  THEN 
PRINT  640,  J 
READ (5,  570)  LINE 
L I NLABEL ( J ) *L I NE 
END  IF 
END  DO 

CALL  WAVORAF(  X.  0.  0.  X  ( L  ) .  AR,  O  ,  1.  ,  L,  NW,  3000,  15,  LINLABEL, 
RXLABEL. YLABEL.  TITLE.  .  TRUE.  .  .  TRUE.  ) 

GO  TO  500 
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PLOT  PARTICLE  DISPLACEMENT  TIME  PROFILE 

IF  < NW.  EQ.  1)  CO  TO  117 
DO  I  =  1 > NW 

LINLABEL(I)-'  ' 

END  DO 
J-l 

DO  WHILE  <  J  .  NE.  O  > 

PRINT*. 'ENTER  THE  LINE  *  TO  BE  LABELED' 

PRINT*, '(TYPE  0  IF  NO  LABELS  NEEDED  >' 

READ  (5.*)  J 
IF  (  J  .  NE.  0  )  THEN 
PRINT  640.  J 
READ <5,  570)  LINE 
L I NLABEL  <  J ) =L I NE 
END  IF 
END  DO 

TITLE* '*P ARTICLE  DISPLACEMENT  PROFILE*' 

XLABEL* '*T ime  (secXl.E6)*' 

YLABEL* '*Disp laccment  (inXl. E6)*' 

DPLOT  ( 1  >=0.  0 

TP ( 1 )=T(NV,  I >*1.  E6 

DO  1=2, NW 

TP ( I )  =  T(NV,  I )*1,  E6 
DPLOT < I ) =DISP ( NV,  I)*i.  E6 
END  DO 
DISPMAX=0.  O 
DO  1=1, NW 

IF  ( DPLOT ( I ) . GE. DISPMAX )  DISPMAX=DPLOT ( I ) 

END  DO 

WR ITE ( 6, 1050)  DISPMAX 

PRINT*. 'ENTER  MAX.  ORDINATE  OF  DISPLACEMENT' 

READ <  5.*)  DHI 
DLO  =00 
TMIN=0.  0 

TMAX=T(NV, NW ) * 1 .  E06 

CALL  WAVGRAF (TP,  TMIN, TMAX. DPLOT,  DLO,  DHI,  NW,  1,  15.  1, LINLABEL. 
ttXLABEL,  YLABEL.  TITLE.  .  TRUE.  .  .  TRUE.  ) 

GO  TO  500 

PLOT  PARTICLE-VELOCITY  TIME  PROFILE 

IF  (NW  EQ  1)  GO  TO  117 
DO  1=1, NW 

LI NLABEL ( I )= '  ' 

END  DO 
J-l 

DO  WHILE  <  J  .  NE  0  ) 

PRINT*. 'ENTER  THE  LINE  #  TO  BE  LABELED' 

PRINT*. '(TYPE  0  IF  NO  LABELS  NEEDED  )' 

READ  (5,*)  J 
IF  (  J  .  NE.  0  )  THEN 
PRINT  640,  J 
READ (5, 570)  LINE 
L I NLABEL ( J  >  =L I NE 
END  IF 
END  DO 

TITLE- '*PARTICLE  VELOCITY  PROFILE*' 

XLABEL- '*Tim»  (secXl  E6>*' 

YLABEL- '*V»locity  (in/s)*' 
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VP  ART  <  1  >=0.  0 

TP ( 1 >=T(NV, 1 )*1E6 

DO  1*2, NW 

TP  ( I )  =  T(NV,  I)*l.  E6 
VP ART ( I )=V(NV, I ) 

END  DO 
VMAX=0.  0 
DO  1*1.  NW 

IF  <  VPART ( I ) .  GE.  VMAX )  VMAX*VPART < I ) 

END  DO 

WRITE (6,  850)  VMAX 

PRINT#, 'ENTER  THE  MAXIMUM  ORDINATE  OF  VELOCITY' 

READ  (5,*)  VHI 
VL0*0.  O 
TMIN*0  O 

TMAX*T  <NV» NW ) * 1 .  E6 

CALL  WAVGRAF(TP,  TMIN,  TMAX,  VPART. VLO, VHI, NW,  1.15,  1.  L INLABEL. 
ttXLABEL,  YLABEL,  TITLE.  .  TRUE.  .  .  TRUE.  ) 

PRINT*,  'INTERESTED  IN  ANOTHER  PROFILE  PLOT.  ENTER  Y  OR  N ' 
READ( 5,  510)  ANS 
IF  (ANS.  EG.  'Y'>  GO  TO  75 
GO  TO  500 

PLOT  PARTICLE  ACCELERATION  PROFILE 

IF  (NW  EQ.  1 )  CO  TO  117 
DO  J*l,  NW— 1 

AB ( 2*J— 1 )=A(NV,  J)/1000. 

AB ( 2* J ) *A ( NV,  Jl/1000. 

LINLABEL(J)*'  ' 

•  END  DO 
J»1 

DO  WHILE  (  J  .  NE  0  ) 

PRINT*.  'ENTER  THE  LINE  #  TO  BE  LABELED' 

PRINT*. '(TYPE  O  IF  NO  LABELS  NEEDED  >' 

READ(  5.  *)  J 
IF  (  J  .  NE.  0  )  THEN 
PRINT  640.  J 
READ (5. 570)  LINE 
L I NLABEL ( J ) *L I NE 
END  IF 
END  DO 

TPLOT ( 1 ) =T ( NV,  1 )*1E6 
DO  J*l. NW— 2 

TPLOT (2*J)-T(NV.  J+l )*1.  E6 
TPLOT ( 2*U+1 >*T(NV.  J+l )*1.  E6 
END  DO 

TPLOT ( 2*NW— 1 ) ■ <  T ( NV.  NW-1 >+DELTR)*l.  E6 
MNW=NW*2 
AMAX*0.  0 
DO  1*1, MNW 

IF  (AB(NV,  I ).  GE.  AMAX)  AMAX=AB (NV,  I ) / 1000. 

END  DO 

WRITE (6, 1100)  AMAX 

PRINT*,  'ENTER  THE  MAXIMUM  ORDINATE  OF  AMPLITUDE' 

READ  (5,*)  AH  I 

TITLE*  '*P ARTICLE  ACCELERATION  PROFILE* ' 

XLABEL* '*Time  (»ecXl.E6)*' 

YLABEL* '*Acc*l»r«t ion  (kin/s/s)*' 

TMIN-0  0 


o  o  o  non 
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TMAX*T(NV, NW)*1.  E6 

CALL  WAVGRAF < TPLOT .  TMIN.  TMAX,  AB.  ALO.  AHI.  MNW.  1. 40,  i. L INLABEL, 
ttXLABEL, YLABEL. TITLE,  .  TRUE.  ,  .  TRUE.  ) 

PRINT*, 'INTERESTED  IN  ANOTHER  PROFILE  PLOT.  ENTER  Y  OR  N' 

READ ( 5, 510)  ANS 
IF  <ANS.  EG.  *Y*>  GO  TO  75 
CO  TO  500 
117  WRITE(6.  550) 

GO  TO  500 

PLOT  THE  STRESS  TIME  PROFILE 

125  IF  <NW. EG. 1)  CO  TO  117 
DO  J= l ,  NW 

LINLABEL( J>-'  ' 

PSTRES ( J  > -STRES  <  NV, J ) 

TP(J)«T(NV>  J>*1.  E6 
END  DO 
STRMAX=0. 0 
DO  I  *»  1 »  NW 

IF  (STRES(NV,  I ).  GE.  STRMAX)  STRMAX-STRES(NV.  I > 

END  DO 

WRITE(6. 950)  STRMAX 
STRL0=0.  O 

PRINT*, 'ENTER  THE  MAXIMUM  ORDINATE  OF  STRESS' 

READ < 5.*)  STRHI 
TITLE* '*STRESS  PROFILE*' 

XLABEL*  '*T  ini*  (secXl.E6)*' 

YLABEL*  '*Str«ss  (psi)*' 

J-l 

DO  WHILE  (  J  .  NE  0  > 

PRINT*. 'ENTER  THE  LINE  *  TO  BE  LABELED' 

PRINT*. '(TYPE  O  IF  NO  LABELS  NEEDED  )' 

READ  <  5.  *)  J 
IF  <  J  .  NE.  0  )  THEN 
PRINT  640.  J 
READ (5. 570)  LINE 
L I NLABEL  ( J ) =L I NE 
END  IF 
END  DO 
TMIN=0  0 

TMAX=T(NV,  NW ) * 1 .  E6 

CALL  WAVGRAF < TP. TMIN. TMAX. PSTRES,  STRLO,  STRHI.  NW,  1.  15.  1, 

«LINLABEL. XLABEL. YLABEL. TITLE.  .  TRUE.  .  .  TRUE.  ) 

PRINT*. 'INTERESTED  IN  ANOTHER  PROFILE  PLOT,  ENTER  Y  OR  N' 

READ! 5. 510)  ANS 

IF  (ANS.  EG.  'Y'>  GO  TO  75 

GO  TO  500 

PROBABILISTIC  CALCULATIONS 
135  IF  (NR.  EG.  1)  THEN 

PRINT*. 'CANNOT  REQUEST  PROBABILISTIC  CALCULATIONS  WITH  NR-1 ' 

GO  TO  500 
END  IF 

190  PRINT*,  '  ' 

PRINT*. 'ENTER  DEPTH  PARAMTER  <  N  >  FOR  THE  PROBABILISTIC  PLOTS' 
PRINT*.  '  ' 

READ  (5.  •>  NV 
XCURNT«NV*QL/40 
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PRINT  680,  XCURNT 
DO  J  -  1 >  NW 

TAVC< J>-( < T 1 < NV,  J)-*-T2(NV. J)*T3(NV,  J)+T4(NV.  J) >/4.  )*1.  E6 
AMPLITUDE  CALCULATIONS 

EA(J)*( A1 (NV,  J ) +A2 < NV*  J ) +A3 < NV,  U)+A4(NV.  J) >/4. 

EA2*(A1(NV.  J>*#2+A2(NV.  J ) *#2+A3 ( NV,  J > ##2+A4 ( NV,  J>#*2>/4. 
SA-SQR T ( EA2-EA <  J ) **2 ) 

AP  ( J )  =E  A  <  J )  +SA 
AM  <  J ) =EA  <  J ) -SA 

STRESS  CALCULATIONS 

ES< J>-(S1(NV, J)+S2(NV,  J)+S3(NV.  J>+S4(NV,  J))/4. 

SPLOT ( J,  1 ) =ES<  J) 

ES2=(S1(NV. J)**2+S2(NV. J>**2+S3(NV,  J ) *#2+S4 ( NV,  J)**2)/4. 
SS-SQRT ( ES2-ES ( J ) **2 ) 

SP( J)*ES( J)+SS 
SPLOT < J,  2)«SP( J> 

SM< J)»ES< J)-SS 
SPLOT (J, 3)*SM( J) 

VELOCITY  CALCULATIONS 

EV( J)*(V1(NV,  J ) +V2 < NV,  J)+V3(NV.  J)+V4(NV.  U) >/4. 

EV2=( VI  (NV,  J)**2+V2(NV.  J)**2-*-V3(NV»  U  >  **2«-V4  (  NV.  J>*#2>/4.  O 
VV< J, 1>=EV( J) 

SV*SQRT  <  EV2-EV <  J ) **2 ) 

VP( J)*EV(J)+SV 
VV(J.  2)=VP(J) 

VM< J)=EV< J)-SV 
VV(J.  3  )  =VM  ( J) 

DISPLACEMENT  CALCULATIONS 

ED( J)*(D1 (NV,  J ) +D2 < NV,  J)+D3(NV,  J)+D4(NV.  J) )/4.  O 
DD< J.  1 >-ED( J>*1.  E6 

ED2=(D1  (NV,  J)**2+D2(NV,  J)**2+D3(NV.  J)**2+D4(NV.  J)**2)/4.  O 
SD=SQRT ( ED2-ED ( J ) **2 ) 

DP(J)*ED( J)+SD 
DD( J.  2)»DP(J)*1.  E6 
DM  <  J ) *ED ( J ) -SD 
DD( J, 3>-DM( J)*l.  E6 
END  DO 

MODIFY  AND  RE-ASSION  THE  AMPLITUDE  AND  TIME  ARRAYS 
DO  1*1,  NW- 1 

APLOT ( 2*1-1 .  1  >-EA(  D/1000. 

APLOT ( 2* 1 ,  1 >=EA< I )/10O0. 

APLOT ( 2*I~ 1 , 2)*AP ( I ) /1000. 

APLOT ( 2*1 , 2>*AP ( I ) / 1000. 

APLOT ( 2*1-1 ,  3 ) "AM ( I )  / 1 000. 

APL0T(2*I. 3>-AM( I )/1000. 

L INLABEL ( I ) * '  ' 

END  DO 

TPL0T(1)»TAV0(1) 

DO  1*1 , NW— 2 

TPLOT (2*1) -TAVO ( I ♦ 1 ) 


TPLOT  <2*1  +  1 ) -TAV© <  I  + 1 ) 

END  DO 
TMIN-O. 0 

TMAX-T  <NV»  NW-1 )*1.  E06 

TPLOT ( 2*NW-1 )-TAVG(NW-l)+DELTR*l.  E6 

ALO-O  O 

DLO-O  0 

VLO-O.  0 

SLO-O.  0 

PRINT*,  '  ' 

PRINT*, 'AVAILABLE  PROBABILISTIC  OPTIONS  ARE' 

PRINT*.  '  ' 

PRINT*. '1.  AMPLITUDE  PROFILE  PLOT' 

PRINT*.  '  ' 

PRINT*. '2.  VELOCITY  PROFILE  PLOT' 

PRINT*. '  ' 

PRINT*. '3.  STRESS  PROFILE  PLOT' 

PRINT*.  '  ' 

PRINT*,  '4.  DISPLACEMENT  PROFILE  PLOT' 

PRINT*. '  ' 

PRINT*. 'ENTER  THE  OPTION  NUMBER' 

READ (5.  580)  JF 
GO  TO  (210.220.230,235)  JF 
210  TITLE-  ♦PROBABILISTIC  AMPLITUDE  PROFILE*' 

XLABEL='*Tim*  (secXl.E6)*' 

YLABEL- '*AMPL I TUDE  <kin/»/s)«' 

APL-O.  O 
NNW»<NW-1)*2 
DO  J-l,  3 
DO  1  =  1,  NNU 

•  IF' (APLOT(I.  J)  .  CE.  APL)  APL-APLOT ( I ,  J > 

END  DO 
END  DO 

PRINT  1110,  APL 

PRINT*.  'ENTER  MAXIMUM  ORDINATE  OF  PROBABILISTIC  AMPLITUDE' 

READ ( 5. *)  AH I 

J=1 

DO  WHILE  (  J  .  NE.  O  ) 

PRINT*, 'ENTER  THE  LINE  #  TO  BE  LABELED' 

PRINT*. '(TYPE  0  IF  NO  LABELS  NEEDED  >' 

READ  (5.  *)  J 
IF  (  J  .  NE.  0  )  THEN 
PRINT  640.  J 
READ (5. 570)  LINE 
L I NLABEL ( J  >  =L I NE 
END  IF 
END  DO 

CALL  WAVGRAF( TPLOT,  TMIN,  TMAX,  APLOT.  ALO.  AHI,  NNW,  3,  30.  3. 

#L INLABEL,  XLABEL,  YLABEL,  TITLE,  .  TRUE.  ,  .  TRUE.  > 

GO  TO  520 
C 

220  TITLE- '^PROBABILISTIC  VELOCITY  PROFILE*' 

XLABEL- '*Ti m*  (stcXl.E6)*' 

YLABEL- ' «V* 1 oc i t g  (in/»»e)*' 

VPL-O.  0 
DO  J— 1 »  3 
DO  I-l.NW 

IF  (VV(I,  J).  OE.  VPL)  VPL-W(I.J) 

END  DO 
END  DO 
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PRINT  1120,  VPL 

PRINT*.  'ENTER  MAXIMUM  ORDINATE  OF  PROBABILISTIC  VELOCITY' 

READ (5,*)  VHI 

J=1 

DO  WHILE  (  J  .  NE.  0  > 

PRINT*, 'ENTER  THE  LINE  tt  TO  BE  LABELED' 

PRINT*, '(TYPE  0  IF  NO  LABELS  NEEDED  )' 

READ <5, *>  J 
IF  (  J  .  NE.  0  )  THEN 
PRINT  640,  J 
READ (5.  570)  LINE 
LINLABEL(J)«LINE 
END  IF 
END  DO 

CALL  WAVORAF ( TAVQ,  TMIN,  TMAX,  VV,  VLO.  VHI,  NW.  3,  15,  3.  LI ML ABEL 
ft, XLABEL, YLABEL,  TITLE.  .  TRUE.  ,  .  TRUE.  > 

CO  TO  520 

C 

230  TITLE® '^PROBABILISTIC  STRESS  PROFILE*' 

XLABEL® '*T iffl*  <s»cXl.E6>«' 

YLABEL® '*Str**s  (ksi)*' 

SPL=0.  0 
DO  J*  1,3 
DO  I ® 1 ,  NW 

IF  (SPLOT(I,  J).  GE.  SPL)  SPL*SPLOT< I,  J> 

END  DO 
END  DO 

PRINT  1130.  SPL 

PRINT*,  'ENTER  THE  MAXIMUM  ORDINATE  OF  STRESS' 

READ  (5.*)  SHI 
J®1 

DO  WHILE  (  J  NE.  0  )  « 

PRINT*, 'ENTER  THE  LINE  #  TO  BE  LABELED' 

PRINT*, '(TYPE  0  IF  NO  LABELS  NEEDED  )' 

READ <  5,  *)  J 
IF  (  J  .  NE.  0  )  THEN 
PRINT  640.  J 

I  READ (5. 570)  LINE 

LINLABEL(J)-LINE 
END  IF 
END  DO 

CALL  WAVGRAF ( TAVG.  TMIN.  TMAX.  SPLOT,  SLO.  SHI, NW, 3,  15, 3.  L INLABEL 
#, XLABEL,  YLABEL.  TITLE,  .  TRUE.  ,  .  TRUE.  ) 

GO  TO  520 

235  TITLE® ^PROBABILISTIC  DISPLACEMENT  PROFILE*' 

XLABEL® '*Tim»  (**cXi.E6)*' 

YLABEL® '*Di*plftCft«*nt  (inXl.  E6)*' 

DPL*0.  O 
DO  J»l,  3 
DO  I ® 1 ,  NW 

IF  (DD(I.  J).  OE.  DPL)  DPL®DD( I , J) 

END  DO 
END  DO 

PRINT  1140,  DPL 

PRINT*.  'ENTER  THE  MAXIMUM  ORDINATE  OF  DISPLACEMENT' 

READ (5,  *)  DHI 
J»1 

DO  WHILE  (  J  .  NE.  0  ) 

PRINT*.  'ENTER  THE  LINE  •  TO  BE  LABELED' 

PRINT*,  '(TYPE  0  IF  NO  LABELS  NEEDED  )' 


A-28 


READ  (5,  *)  J 
IF  <  J  .  NE.  O  )  THEN 
PRINT  640.  J 
READ 0,570)  LINE 
L I NLABEL  ( J ) *L I NE 
END  IF 
END  DO 

CALL  WAVCRAF < TAVC.  THIN.  THAX.  DD.  DLO.  DHI. NW. 3.  15. 3.  LINLABEL 
#. XLABEL, YLABEL. TITLE,  .  TRUE.  .  TRUE.  ) 

GO  TO  520 

PLOT  THE  VOLUME  DISTRIBUTION  FUNCTION 

0  TITLE* '*VOLUME  DISTRIBUTION  FUNCTION*' 

XLABEL='*Dapth  <in>*' 

YLABEL*  '*No  Function*' 

J-l 

DO  WHILE  (  J  .  NE.  O  > 

PRINT*. 'ENTER  THE  LINE  «  TO  BE  LABELED' 

PRINT*.  '(TYPE  0  IF  NO  LABELS  NEEDED)' 

READO.  *)  J 
IF  <  J  .  NE.  O  )  THEN 
PRINT  640.  J 
READO,  570)  LINE 
L I NLABEL  ( U  >  *L I NE 
END  IF 
END  DO 

CALL  WAVGRAF(X,  O.  0.  X(L).  C3.  .  5.  1.  O,  L,  1.  3000.  15,  LINLABEL, 
»XLABEL,  YLABEL.  TITLE.  .  TRUE.  ,  .  TRUE.  ) 

GO  TO  500 
•  * 

PLOT  THE  MU  COEFFICIENT  OF  AMPLITUDE 

140  C 1MIN=200. 

ClMAX=-200. 

DO  J-l.NW 
DO  1*1. L 

IF  (Cl(I.J)  LE.C1MIN)  C1MIN*C1  <  I.  J) 

IF  (Cl < I. J) .  GE.  Cl MAX )  C1MAX*C1(I.  J) 

END  DO 
END  DO 

WRITE(6. 700)  C1MIN, C1MAX 

PRINT*,  ENTER  LOWER  &  UPPER  Y-COORDINATES  OF  THE  GRAPH' 
READO.*)  YLO.  YHI 

TITLE* '*MU-COEFFICIENT  OF  AMPLITUDE*' 

XLABEL* '*Propagation  Distance  (in)*' 

YLABEL* '*Mu-Co»f f iciant  (l/in)*' 

CALL  WAVGRAF ( X.  0  0.  X(L).  Cl.  YLO,  YHI, L,  1. 3000,  15, LINLABEL, 
4XLABEL. YLABEL. TITLE.  .  TRUE.  ,  .  TRUE.  ) 

GO  TO  500 

PLOT  THE  KAPPA  COEFFICIENT  OF  AMPLITUDE 

150  C2M IN-200 
C2MAX— 200 
DO  J-l.NW 
DO  1*1, L 

IF  (C2( I.  J) .  LE.  C2MIN)  C2MIN»C2( I.  J) 

IF  (C2(I, J).  GE.  C2MAX)  C2MAX-C2( I,  J) 

END  DO 


o  o  n  n  o 


END  DO 

WRITE(6. 750)  C2MIN. C2MAX 

PRINT*.  'ENTER  LOWER  &  UPPER  Y-COORDINATES  OF  GRAPH' 

READ (5. *)  YLO.  YHI 

TITLE*' *KAPP A-COEFF I C I ENT  OF  AMPLITUDE*' 

XLABEL* '*Prop«gation  Distanc*  <ln)*' 

YLABEL* '*K— Coaf f ic iant  (sac  2/in  2)*' 

CALL  WAVGRAF ( X > O.  O.  X(L).  C2.  YLO.  YHI.  L.  1. 3000.  19.  L INLABEL. 

4XLABEL,  YLABEL,  TITLE,  .  TRUE.  .  .  TRUE.  > 

GO  TO  500 

MAIN  ROUTINE  COMPUTATION  AND  GRAPHICS  ENDS  HERE 

WRITE  OUTPUT  ON  A  DATA  FILE 

200  OPEN  <  UN I T- 1 .  F I LE= ' DATA ' .  STATUS* ' NEW ' ) 

WRITEd.  240) 

240  FORMAT(  IX.  'N',  7X,  'X(N)  '.  7X.  'UREAL',  7X,  'U(N)  ',  7X,  'A'.  7X.  'AR'. 

#7X,  'Cl  (N,  J) '.  7X.  'C2(N,  J)  '»  7X,  'C3(N,  J)'> 

DO  400  J*l,  NW 
DO  300  N=l,  L 

WRITE(  1,250)  N.  X<N) ,  UREAL  (N,  J),  U<N,  J).  A(N,  J ) ,  AR  ( N,  J),  CUN,  J), 
#C2(N, J),  C3(N,  J) 

250  FORMAT (IX,  14, 8E13.  5/ ) 

300  CONTINUE 
400  CONTINUE 

CLOSE ( UNIT*1 ) 

500  WRITE (6,  *)  'WOULD  LIKE  TO  CONTINUE,  ENTER  Y  OR  N' 

READ (5. 510)  ANS 
510  FORMAT (A3) 

IF  (ANS.  EQ.  'Y'>  GO  TO  72 
IF ( ANS.  EQ.  'N')  GO  TO  630 

IF ((ANS.  NE.  'Y ' > .  AND.  (ANS.  NE.  'N')>  GOTO  500 
520  PRINT*. 'WOULD  LIKE  ANOTHER  PROBABILISTIC  PLOT' 

READ (5. 540)  PANS 
IF(PANS.  EQ.  'Y')  GO  TO  190 
IF(PANS.  EQ.  'N'>  GO  TO  500 

IFdPANS.  NE.  'Y')  AND.  (PANS.  NE.  'N')>  GOTO  520 
540  FORMAT (A3) 

550  FORMAT (IX,  'MORE  THAN  2  AMPLITUDES  NEEDED  TO  MAKE  THIS  PLOY') 

560  FORMAT (21 5) 

570  FORMAT (A20) 

580  FORMAT (112) 

600  FORMAT (IX,  'NOTE:  MAX.  ALLOWABLE  N  FOR  THIS  RUN  IS  -  ',14) 

630  STOP 

640  FORMAT (IX,  'ENTER  THE  LABEL (* — •)  OF  LINE',  12) 

650  FORMAT (115) 

680  FORMAT (IX,  'CURRENT  DEPTH  IS  X  *',FB.  5,  'inch') 

700  FORMAT (IX,  'MINIMUM  MU  *',F10.  5/ 

#1X,  'MAXIMUM  MU  *',F10.  5) 

750  FORMAT (IX,  'MINIMUM  KAPA  *',F10.  5/ 

«1X.  'MAXIMUM  KAPA  *',F10.  5) 

760  FORMAT (IX,  'ENTER  THE  <N>  VALUE  CORRESPONDING  TO  DEPTH  REQUIRED') 
800  FORMAT (IX,  'AMPLITUDE( '<  12,  ')  *',E13.5> 

050  FORMATdX.  'MAX.  PARTICLE  VELOCITY  FIO.  5) 

900  FORMATdX.  'MINIMUM  ACTUAL  SPEED  *',F10.  5/ 

#1X.  'MAXIMUM  ACTUAL  SPEED  *',F10.  9) 

950  FORMATdX.  'MAXIMUM  8TRES3  *',F10.  4) 

1000  FORMATdX,  'MAXIMUM  AVERAOE  SPEED  *'.F10.  5,/> 

1050  FORMATdX.  'MAXIMUM  PARTICLE  DISPLACEMENT  *',FB.  4/> 


non  n  n  n  non  non  non  r>  o  o  o  o 


1100  FORMAT <  1 X.  'MAXIMUM  DETERMINISTIC  AMPLITUDE  *',F10.  5, /> 
1110  FORMAT <  1 X.  'MAXIMUM  PROBABILISTIC  AMPLITUDE  -',F8.  4/> 
1120  FORMAT (IX.  'MAXIMUM  PROBABILISTIC  VELOCITY  -'.FB.  4/> 

1130  FORMAT (IX.  'MAXIMUM  PROBABILISTIC  STRESS  -'.FB.  4/) 

1140  FORMAT (IX.  'MAXIMUM  PROBABILISTIC  DISPLACEMENT  -'.FB.  4/) 
END 

SUBROUTINE  GAUSS  (  ENUA.  JF.  KK.  NF.  ZL.  X.  UBAR  ) 


GAUSS-LEGENDRE  QUADRATURE  FOR  CALCULATION  OF  PROPAGATION  TIME  • 
AND  WAVE  VELOCITY  USING  FOUR -POINT  FORMULA  FOR  INTEGRATION  * 


REAL  ETA, FKCH. ENUA, ENUB. NUZX, UF. FXCH, ENU. ENUX.  X 
DIMENSION  ETA<4), NUZX < 4), FXCH< 4 ) , ENU< 4 ) . ENUX ( 4)  . 

#UF ( 4 ) ,  FKCH<  4 ) ,  PENU<4) , PNU<4> 

COMMON  /GRAIN/  AFFO, AFFF,  BB,  DUM, ENUB, FLFFO,  FLFFF,  GO.  GL,  GO. 
#TIME 

DEFINE  THE  CONSTANTS  AND  THE  WEIGHTING  COEFFICIENTS 

PI-3  141592741012573 
WT1-0.  34785485 
WT2-0  65214515 
ETA<  1 )—  O  86113631 
ETA<2>=-0  33998104 
ETA<4)— ETA<  1  ) 

ETA<3>—  ETA<2) 

DEFINE  THE  FUNCTION  AND  PERFORM  THE  INTEGRATION 
DO  60  1  =  1.4 

FKCH ( I >=<ETA< I )*<X-ZL)  +  <X+ZL) )/2.  O 
CO  TO  <20. 30. 40)  KK 

PERIODIC  VOLUME  DISTRIBUTION  FUNCTION 

20  ENU  < I ) -ENUA+  < 1 .  -ENUA ) *COS  <  2.  *P I *FKCH < I > /GL ) 

NUZ X ( I )  —  <  < ENUA- 1 .  0)/GL)*2.  0*PI*SIN<2.  #PI#FKCH< I >/GL> 
PP-FLFFO+O.  5*AFF0*NUZX  < I ) **2 
IF<PP  LT.  0  O)  GO  TO  150 
UF ( I ) -SORT  <  PP ) 

FXCH( I )=1. O/UF ( I > 

GO  TO  50 

EXPONENTIAL  VOLUME  DISTRIBUTION  FUNCTION 

30  PENU < I >  =1 .  0-<l  O-ENUB ) *EXP ( -BB*CG*FKCH< I ) ) 

ENUX  < I ) - < 1  O-ENUB ) *BB*GG*EXP  < -BB*GC*FKCH< I ) ) 

TEST -FLFFO+O.  5*AFF0*ENUX < I > *#2 

IF < TEST  LT  0)  GO  TO  150 

UF < I ) -SORT  < FLFFO+O  5*AFF0*ENUX  < I ) **2 ) 

FXCH< I )  =  1.  O/UF ( I ) 

GO  TO  50 

PERIODIC  AND  EXPONENTIAL  COMBINED  VOLUME  DISTRIBUTION  FUNCTION 

40  ENU( I ) -ENUA+  < 1 .  -ENUA)*C0S<2.  *PI*FKCH< I >/GL) 

NUZX ( I )  =  ( (ENUA-l .  0 ) /GL ) *2.  0*PI*SIN(2.  *PI*FKCH< I >/0L> 

PENU< I )-l  0-< 1  O-ENUB >*EXP < -BB*00*FKCH< I > ) 

ENUX  <!)-<!  O-ENUB )*BB*0G*EXP(-BB*00*FKCH( I ) ) 


nnn  nnn  nnn  nnnnn  n  n  n  n  n  n  n  r> 


PNU< I ) -ENU ( I >*ENUX ( I )+PE NU< I ) *NUZX  < I > 

CHEK-FLFFO+.  S*AFFO*PNU< 1 )**2 

IFtCHEK.  LT.  0>  00  TO  ISO 

UF  <  I )  *  SORT  <  FLFF0+0.  5*AFF0*PNU  < I ) **2  > 

FXCH<I)-1.  0/UF < I ) 

50  CONTINUE 

100  TAUX  =  WT1*(FXCH< 1 )+FXCH(4) )+WT2*<FXCH(2)+FXCH<3) > 

TIME  *  TAUX*( <X-ZL)/2.  0) 

TIME  =  DUM  +  TIME 

CALCULATE  THE  AVERAGE  WAVE  SPEED 

UBAR  »  X/TIME 
GOTO  200 

150  WRITE<&»  180)  NF.  UF 

160  FORMAT ( IX*  ' INPUT  IS  INCOMPATIBLE  WITH  REAL  WAVE  SPEED  IN  GAUSS'/ 
#1X,  'LOOPING  PARAMETER  N  -',15. 

#1X,  'WAVE  AMPLITUDE  NUMBER  J  -',13) 

STOP 

200  RETURN 
END 

SUBROUTINE  RUNKUT< AMP, ENUA.  H,  K,  N,  RKAPA,  RMU,  RVDFUN.  RWS.  X.  VO) 


SUBROUTINE  RUNGE-KUTTA  CALCULATES  THE  WAVE  AMPLITUDE  • 

IN  THE  GRANULAR  MEDIA  * 


REAL  AMP, BB, ENU, ENUX.  ENUXX,  ENUZ,  ENUZZ,  H,  PENU,  RKAPA.  RMU,  RVDFUN, 
«RWS,  X,  YO 

COMMON  /GRAIN/  AFFO, AFFF,  BB.  DUM.  ENUB,  FLFFO,  FLFFF,  GO.  OL.  GO  TIME 
COMMON  /CONTRL/  P.Pl.Q.  01 
DATA  PI/3.  141592741012573/ 

DEFINE  THE  VOLUME  DISTRIBUTION  FUNCTIONS  AND  THEIR  DERIVATIVES 

PERIODIC  VOLUME  DISTRIBUTION  FUNCTION 

ENU <  X  )  —  ( ENUA+  <  1 .  —ENUA )  *COS ( 2.  *PI*X/OL)  )*P«-P1 

ENUZ  <  X )  — <-2.  *PI*<1.  -ENUA)*SIN(2.  *PI*X/OL > /CL > *P 

ENUZ Z ( X ) * ( —4.  *<PI**2.  )*<1.  -ENUA > *COS ( 2.  *PI*X/GL)/<0L**2) )*P 

EXPONENTIAL  VOLUME  DISTRIBUTION  FUNCTION 

PENU(X)*< 1.  0-< 1.  O-ENUB  >*EXP ( -BB*OG*X  > ) *0+01 

ENUX  <  X )  =  < ( 1 .  O-ENUB ) *BB*CO*EXP < -BB«OG*X ) >  *0 

ENUX X  <  X ) - ( - < 1 .  O-ENUB ) *BB**2*GC**2#EXP ( -BB*OQ*X ) ) #Q 

PERIODIC  AND  EXPONENTIAL  VOLUME  DISTRIBUTION  FUNCTION  COMBINED 

PNU( X ) -ENU  <  X ) #PENU( X ) 

PNUX  <  X  >-ENU<  X  )*ENUX  (  X  >-*-PENU<  X  )  *ENUZ  <  X ) 

PNUXX ( X ) =ENU< X  > *ENUXX ( X ) ♦PENU ( X ) *ENUZZ (X)+2  *ENUZ ( X ) *ENUX ( X ) 
COMPUTE  THE  WAVE  AMPLITUDE  AND  THE  CONSTANTS  MU  S<  KAPA 
UFSQ <  X ) - <  FLFFO+O.  5*AFF0*PNUX  <  X ) *#2 ) 

FMU  <  X )  —  <  UFSQ <  X ) /PNU  <  X ) +AFFO*PNUX  X ( X  > /2.  0) *PNUX ( X ) / <2.  *UFSO(X> ) 
CKAP < X >—— (FLFFF ♦AFFF *PNUX < X ) **2/2.  >/<2.  *00*UFS0(X)#*2) 

FC  <  X. Z )— CKAP  <  X ) *Z**2-FMU<  X  >  *Z 
AK1-FC(X,  YO) 
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<J  o 


AK2-FC ( X+H/2.  0.  Y0+H*AKl/2.  0) 

AK3-FC ( X+H/2.  0-  YO+H*AK2/2.  0) 

AK4-FC  <  X +H,  Y0+H*AK3 ) 

Y0  «  Y0+H*(AKl+2.  0*AK2+2.  0*AK3+AK4) /6.  0 

AMP-YO 

RMU-FMU(X) 

RKAPA-CKAP < X > 

RVDFUN-PNU(X) 

IF  <UFSQ(X>.  LT.  0)  CO  TO  10 
R WS-SQR T ( UFSQ ( X ) ) / 1 000 
00  TO  30 

10  WRITE<6.  20)  N,  K 

20  FORMAT (IX, 'INPUT  IS  INCOMPATIBLE  WITH  REAL  WAVE  SPEED  IN  RUNKUT'/ 
#1X,  'LOOPING  PARAMETER  N  -'.15, 

#1X.  'WAVE  AMPLITUDE  NUMBER  J  -',I3> 

STOP 

30  RETURN 
END 


SUBROUTINE  WAVCRAF(X, XMININP.  XMAXINP,  Y,  YMININP,  YMAXINP,  N.  M. 

-  NORIG.  MORIG,  LINLABEL.  INDAX,  DEPAX.  TITLE.  INDEVEN. 

-  DfePEVEN) 

C  ************************************************************* 

c 

C  To  draw  a  lina  graph  with  muitipla  plots  on  an  BX11  paga. 

C  using  DI-3000  graphics  subroutinas 

C 

C  N:  Number  of  values  in  each  set 

C  M:  Number  of  sets  of  dependant  values 

C  NORIG  Original  column  dimension  of  the  value  arrays 

C  MORIG  "  row  M  . 

C  LINLABEL.  Character  string  labels  of  each  ‘.et  of  values 
C  INDAX  Character  string  of  independent  axis  label 

C  DEPAX:  ••  "  -  dependent 

C  TITLE:  "  "  title  of  graph 

C  (Note:  all  character  strings  should  be  inclosed  in  a 

C  delimiter  mark,  such  as  '4graph*') 

C  INDEVEN:  Logical  variable  which  is  .TRUE,  if  the  independent 
C  axis  should  be  extended  to  terminate  at  even  divisions. 

C  DEPEVEN:  Logical  variable  which  is  .TRUE,  if  the  dependent 
C  axis  should  be  extended  to  terminate  at  even  divisions. 

C  VECTOR:  Array  for  graph  information  storage 
C  VSI2E:  Si x e  of  VECTOR  array 
C  RATIO  Aspect  ratio  of  the  graphics  device 
C  NSET:  Data  set  counter 

C  XDIVMAX:  Maximum  number  of  independent  axis  divisions 

C  YDIVMAX :  Maximum  number  of  dependent  axis  divisions 

C  SPAN:  Range  of  axes  (maximum  value  minus  minimum  value) 

C  ORDER:  order  of  the  axis  increments 

C  RESOLUTION:  Minimum  even  incremented  resolution  of  the  axis 

C 

REAL  X<  NORIC  ).  XMININP,  XMAXINP,  Y<  NORIO,  MORIG  ).  YMININP. 

-  YMAXINP 

INTEQER  N,  M,  NORIG,  MORIO 

CHARACTER  LINLABEL! 10)*25,  INDAX*50. DEPAX*50.  TITLE*50 
LOO I CAL  INDEVEN,  DEPEVEN 

REAL  ENTER,  OUTER,  XMIN.  XMAX,  YMIN.  YMAX,  RATIO,  SPAN,  ORDER. 

•  RESOLUTION.  X INCREMENT.  Y INCREMENT 
INTEGER  VECTOR ( l 5000 ) .  VSIZE,  NSET,  XDIVMAX.  YDIVMAX 
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DATA  VSIZE/15000/.  XDIVMAX/5/.  YDIVMAX/10/ 

XMIN-XMININP 

XMAX-XMAXINP 

YMIN-YMININP 

YMAX*YMAXINP 

Initialize  CRAFMAKER  and  DI-3000 


CALL 

JCHINI < 

.  TRUE.  . 

1  ) 

CALL 

JCHART ( 

VECTOR, 

VSIZE 

> 

CALL 

JASPEK ( 

1,  RATIO  ) 

CALL 

UVSPAC  < 

-1.0,  1. 

0,  -RATIO,  RATIO  > 

CALL 

JGRAPH< 

VECTOR, 

VSIZE, 

1  ) 

!ine 

title  and 

1  text  characteristics. 

CALL 

JXTEXT ( 

VECTOR, 

VSIZE, 

1,  S, 

0.0,  1. 

0,  0  ) 

CALL 

J TXBOX< 

VECTOR, 

VSIZE, 

0,  0, 

1  ) 

CALL 

JSTNOT ( 

VECTOR. 

VSIZE, 

1,  1, 

TITLE  ) 

CALL 

JPONOT  < 

VECTOR. 

VSIZE, 

1.  1, 

500.  0, 

980  0  ) 

CALL 

JTXHCT  < 

VECTOR. 

VSIZE, 

20.  0, 

0.  0,  0. 

0  ) 

If  desired,  redefine  the  dependent  axis  minimum  and 
maximum  for  even  axes 

SPAN  =  YMAX  -  YhIN 
IF  (SPAN.  LT.  IE-6)  THEN 

ORDER  *  LOG 1 0 (  SPAN*1E12/ 12.  0  )  -  L0G10<  FLOAT (  YDIVMAX  )  ) 
ELSE 

ORDER  -  L0G10<  SPAN  )  -  L0Q10<  FLOAT (  YDIVMAX  )  > 

END  IF 

IF  (  ORDER  .  LT.  0.0  )  ORDER  -  ORDER  -1.0 
RESOLUTION  -  100  **  INT<  ORDER  ) 

IF  (  SPAN  /  FLOAT <  YDIVMAX  )  .  LE.  RESOLUTION  )  THEN 
Y INCREMENT  -  RESOLUTION 

ELSE  IF  <  SPAN  /  FLOAT <  YDIVMAX  )  .  LE.  2.  O  *  RESOLUTION  )  THEN 
YINCREMENT  *  2.  0  *  RESOLUTION 

ELSE  IF  (  SPAN  /  FLOAT <  YDIVMAX  )  .  LE.  5.  O  *  RESOLUTION  )  THEN 
YINCREMENT  -  5.  0  *  RESOLUTION 

ELSE  IF  <  SPAN  /  FLOAT <  YDIVMAX  )  .  LE.  10.0  *  RESOLUTION  )  THEN 
YINCREMENT  »  10.  0  *  RESOLUTION 
ELSE 

YINCREMENT  -  20.  0  *  RESOLUTION 
END  IF 

IF  <  INDEVEN  )  THEN 

IF  (  YMIN  .  OE.  0  )  THEN 

YMIN  -  YINCREMENT  *  FLOAT <  INT<  YMIN  /  YINCREMENT  ♦  O.  01  >  ) 
ELSE 

YMIN  -  YINCREMENT  *  FLOAT <  INT<  YMIN  /  YINCREMENT  -  0.  99  )  > 
END  IF 

IF  (  YMAX  .  OE  0  )  THEN 

YMAX  -  YINCREMENT  *  FLOAT <  INT<  YMAX  /  YINCREMENT  *  0.  99  )  ) 
ELSE 

YMAX  -  YINCREMENT  *  FLOAT <  INT<  YMAX  /  YINCREMENT  -0.01  )  ) 
END  IF 
END  IF 

CALL  JBTVAX  (VECTOR,  VSIZE.  1,  1,  YMIN,  YMAX,  DEPAX) 

IF  <  .  NOT.  INDEVEN  )  THEN 
YMIN  -  YMIN  ♦  YINCREMENT 
IF  <  YMIN  .  OE.  0  )  THEN 


YMIN  »  YINCREMENT  *  FLOAT <  INT(  YMIN  /  YINCREMENT  +  0.01  >  ) 
ELSE 

YMIN  ■  YINCREMENT  #  FLOAT (  INT(  YMIN  /  YINCREMENT  -  0.  99  >  > 
END  IF 

YMAX  «  YMAX  -  YINCREMENT 
IF  (  YMAX  .  GE.  0  )  THEN 

YMAX  -  YINCREMENT  *  FLOAT <  INT<  YMAX  /  YINCREMENT  +  0.  99  )  ) 
ELSE 

YMAX  -  YINCREMENT  *  FLOAT (  INT<  YMAX  /  YINCREMENT  -0.01  >  ) 
END  IF 
END  IF 

CALL  JTIC  (VECTOR.  VSIZE.  1.  1.  1.  YMIN,  YMAX,  YINCREMENT) 

CALL  JTCATR  (VECTOR,  VSIZE,  1,  1.  1,  0.  0,  6.  0,  0) 

If  desired,  redefine  the  independent  exit  minimum  and 
maximum  for  even  divisions 

SPAN  =  XMAX  -  XMIN 
IF  (SPAN.  LT.  IE-6)  THEN 

ORDER  »  L0G10 <  SPAN+1E12/12  0  )  -  L0G10(  FLOAT (  XDIVMAX  )  ) 
ELSE 

ORDER  *  LOG 1 0 (  SPAN  )  -  L0C10(  FLOAT (  XDIVMAX  )  ) 

END' IF 

IF  (  ORDER  .  LT  0.0  )  ORDER  -  ORDER  -  1.0 
RESOLUTION  =  10.  0  **  INT(  ORDER  ) 

IF  (  SPAN  /  FLOAT (  XDIVMAX  )  .  LE.  RESOLUTION  )  THEN 
X INCREMENT  «  RESOLUTION 

ELSE  IF  (  SPAN  /  FLOAT (  XDIVMAX  )  .  LE.  2.  O  *  RESOLUTION  )  THEN 
X INCREMENT  »  2.  0  *  RESOLUTION 

ELSE  IF  (  SPAN  /  FLOAT (  XDIVMAX  >  LE.  5.  O  *  RESOLUTION  )  THEN 
.  X INCREMENT  «  5  0  *  RESOLUTION 

ELSE  IF  (  SPAN  /  FLOAT (  XDIVMAX  )  .  LE.  1*  O  *  RESOLUTION  )  THEN 
XINCREMENT  -  10.  0  *  RESOLUTION 
ELSE 

XINCREMENT  *  20.  0  *  RESOLUTION 
END  IF 

IF  (  INDEVEN  >  THEN 

IF  (  XMIN  .  CE.  O  )  THEN 

XMIN  =  XINCREMENT  *  FLOAT (  INT(  XMIN  /  XINCREMENT  +0.01  )  ) 
ELSE 

XMIN  =  XINCREMENT  *  FLOAT (  INT(  XMIN  /  XINCREMENT  -  O.  99  )  ) 
END  IF 

IF  (  XMAX  .  GE.  0  )  THEN 

XMAX  =  XINCREMENT  *  FLOAT (  INT(  XMAX  /  XINCREMENT  +  0.  99  )  ) 

ELSE 

XMAX  =  XINCREMENT  *  FLOAT (  INT(  XMAX  /  XINCREMENT  -0.01  )  ) 
END  IF 
END  IF 

CALL  JSTHAX (  VECTOR,  VSIZE,  1,  2,  XMIN,  XMAX,  INDAX  ) 

IF  (  NOT.  INDEVEN  )  THEN 
XMIN  =  XMIN  +  XINCREMENT 
IF  (  XMIN  .  GE.  0  )  THEN 

XMIN  -  XINCREMENT  *  FLOAT (  INT(  XMIN  /  XINCREMENT  +0.01  )  ) 
ELSE 

XMIN  -  XINCREMENT  *  FLOAT (  INT(  XMIN  /  XINCREMENT  -  0.  99  )  ) 
END  IF 

XMAX  *  XMAX  -  XINCREMENT 
IF  (  XMAX  .  GE  0  )  THEN 

XMAX  -  XINCREMENT  *  FLOAT  (INT  (XMAX  /  XINCREMENT  ♦  0.99)) 


non  ooo  ooo  ooo  non  o  n  r>  ooo  non 


XMAX  -  XINCREMENT  *  FLOAT  ( INT  ( XMAX  /  XINCREMENT  -  0.  01 ) ) 
END  IF 
END  IF 

DESCRIBE  TWO  TICK  MARK  CROUPS  TO  BE  USED  ON  THE  X-AXIS 

CALL  JTIC  (VECTOR,  VSIZE,  1.  2.  2,  XMIN,  XMAX.  XINCREMENT) 

CALL  JTCATR (VECTOR.  VSIZE,  1. 2.  2,  0.  0.  IS.  0,  0) 

ENTER-XMIN+X INCREMENT/2.  O 
OUTER*XMAX-X INCREMENT/2.  0 

CALL  JTIC  (VECTOR,  VSIZE,  1,  2,  3.  ENTER.  OUTER,  XINCREMENT) 
CALL  JTCATR  (VECTOR,  VSIZE,  1.  2,  3,  0.  0,  7.  5.  O) 

CALL  JTCPAT  (VECTOR, VSIZE,  1 , 2, 3,  10,  20) 

□PEN  UP  A  LEGEND 


CALL  JTXBOX (  VECTOR, 
CALL  JSTLCD(  VECTOR, 
CALL  JLGPOS(  VECTOR. 
CALL  JTXHGT (  VECTOR, 


VSIZE,  0.  0,  1  ) 

VSIZE,  1.  '•*'  ) 

VSIZE,  1,  500.0,  860.  0  > 

VSIZE,  17.  0.  0.  0.  0.  O  ) 


PASS  THE  DATA  SETS  TO  GRAFMAKER 


CALL  JRDATA  (VECTOR.  VSIZE,  1.  X,  N) 

CALL  JINDEP  (VECTOR.  VSIZE,  1.  1) 

DO  NSET  -  1.  M 

CALL  JRDATA  (  VECTOR,  VSIZE,  NSET  ♦  1,  Y(  1,  NSET  ) ,  N  ) 
CALL  JDEPEN(  VECTOR.  VSIZE,  1,  NSET.  NSET  +  1  ) 

DEFINE  PLOT  LINE  CHARACTERISTICS 

CALL  JXLINE (  VECTOR,  VSIZE.  NSET.  NSET,  16383,  (NSET-1)*10, 

-  16383  ) 

CALL  JDTATR (  VECTOR.  VSIZE,  1.  NSET.  0,  0,  NSET  ) 

MAKE  A  LEOEND  ENTRY 


IF  (  L INLABEL  <  NSET  )(1:1)  EG.  AND.  LINLABEL(  NSET  )(2:2> 

. NE  )  CALL  JSDLGD (  VECTOR,  VSIZE,  1,  NSET. 

L 1 NLABEL (  NSET  )  ) 

END  DO 


SHOW  CHART 


CALL  JCHSHW  (VECTOR,  VSI ZE.  -0.70,  0.70,  -0.55,  0.55) 

PAUSE  FOR  VIEWING 

CALL  JPAUSE(  1  ) 

Terminate  GRAFMAKER 

CALL  JCHTRM(  .  TRUE.  > 

RETURN 

END 
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